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Introduction 


Approximately  2.5  million  Americans  live  with  epilepsy  and  epilepsy-related  deficits  today,  more 
than  disabled  by  Parkinson  disease  or  brain  tumors.  The  impact  of  epilepsy  in  the  US  is  significant 
with  a  total  cost  to  the  nation  for  seizures  and  epilepsy  of  approximately  $12.5  billion.  Epilepsy 
consists  of  more  than  40  clinical  syndromes  affecting  40  million  people  worldwide.  Approximately 
25  percent  of  individuals  receiving  antiepileptic  medication  have  inadequate  seizure  control; 
however,  80%  individuals  with  medication  resistant  epilepsy  might  be  cured  through  surgery  if  one 
were  able  to  precisely  localize  the  seizure  focus.  The  proposed  research  will  significantly  advance 
our  ability  to  localize  such  foci,  and  thereby  offer  curative  epilepsy  surgery  for  this  devastating 
disease.  Photoacoustic  tomography  (PAT)  uniquely  combines  the  high  contrast  advantage  of 
optical  imaging  and  the  high  resolution  advantage  of  ultrasound  imaging  in  a  single  modality.  In 
addition  to  high  resolution  structural  information,  the  proposed  PAT  is  also  able  to  provide 
functional  information  that  are  strongly  correlated  with  regional  or  focal  seizure  activity,  including 
blood  volume  and  blood  oxygenation  because  of  the  high  sensitivity  of  optical  contrast  to 
oxyhemoglobin  and  deoxyhemoglobin  concentrations.  The  hypothesis  of  the  proposed  research  is 
that  PAT  offers  the  possibility  to  non-invasively  track  dynamical  changes  during  seizure 
occurrence.  The  overall  goal  of  this  research  is  to  advance  a  finite  element  based  photoacoustic 
tomography  method  for  epilepsy  imaging,  using  both  laboratory  and  in  vivo  experiments. 
Specifically,  in  this  project  we  propose:  (1)  To  design,  construct  and  test  a  transducer  array 
system  for  both  2D  and  3D  PAT  imaging;  (2)  To  advance  reconstruction  algorithms  and  associated 
image  enhancement  schemes  for  quantitative  PAT;  (3)  To  evaluate  and  optimize  the  integrated 
functioning  of  the  hardware  and  software  components  of  the  transducer  array-based  system,  using 
simulation  and  phantom  experiments;  (4)  To  test  and  validate  the  PAT  system  using  a  well 
established  animal  model  of  temporal  lobe  epilepsy. 

Body 

This  report  describes  work  accomplished  in  4  years  (plus  1  year  no-cost-extension)  of  the  project. 
As  outlined  in  the  approved  Statement  of  Work  (SOW),  the  tasks  during  this  period  of  time  include: 
Task  1  (Months  0-24):  Purchase  and  calibrate  Ti:Sapphire  tunable  laser;  design  and  test  PVDF 
transducers;  design  and  test  data  acquisition  subsystem;  assemble  the  entire  photoacoustic 
tomography  (PAT)  system;  Task  2.  (Months  18-24)  Design,  build  and  test  the  animal  interface 
Task  3  (Months  0-24):  Implement  reconstruction  codes  and  associated  enhancing  schemes 
including  dynamic  dual  meshing,  total  variation,  advanced  regularization  techniques,  adaptive 
meshing,  initial  parameter  optimization  and  reconstruction  of  absolute  optical  absorption  coefficient 
for  quantitative  high  resolution  functional  PAT.  Task  4.  (Months  12-24)  Calibrate  the  imaging 
system;  conduct  simulation  and  extensive  phantom  experiments  using  the  system  for  evaluation  of 
the  reconstruction  codes  and  associated  enhancements.  Task  5.  (Months  25-42)  Implement 
reconstruction  codes  in  the  areas  of  advanced  boundary  conditions,  3D  algorithms  and 
parallelization  of  3D  codes;  perform  simulation  and  phantom  experiments  to  evaluate  the  software 
developments.  Task  6.  (Months  25-42)  Conduct  animal  experiments  for  evaluating  the  PAT 
imaging  system.  Task  7.  (Months  43-48)  Analysis  of  the  images  from  the  in  vivo  experiments. 

The  sections  below  consist  of  (1)  Hardware  development  (2)  System  Development  and  evaluation 
(3)  Software  Development  (4)  Animal  experiments  (5)  Rat  Model  of  Temporal  Lobe  Epilepsy  (6) 
Analysis  of  the  images  from  the  in  vivo  experiments  that  reflect  the  tasks  associated  with  the  SOW 
during  Months  0-48. 
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1.  Hardware  Development  (Task  1) 

We  have  completed  the  design  of  the  proposed  array  based  PAT  system  based  on  the  extensive 
modeling  and  testing  of  the  key  components  needed  for  the  proposed  system  using  the  existing 
single  transducer  scanning  system. 

In  Fig.  1a,  the  designing  work  for  an  array  based  PAT  system  is  shown. 
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Fig.  1  (a)  Schematic  of  the  PAT  system,  (b)  antj  (c):  3D  (display  of  the 
system/measurement  chamber  with  (b)  and  without  (c)  the  top  coverecd. 

Accomplished  goals  in  the  first  year  of  the  project  are  summarized  as  follows: 

1 .  We  completed  the  design  for  an  array  based  PAT  system,  which  is  the  most  important  part  of 
Task  1 . 

2.  We  developed  three  novel  schemes  that  can  enhance  our  current  reconstruction  software. 

3.  We  conducted  phantom  experiments  that  confirmed  our  software  enhancement. 

For  more  details  please  see  the  full  progress  report  of  the  first  year. 

2.  System  Development  and  Evaluation  (Tasks  1,  2,  3  and  4) 

In  the  second  year,  we  achieved  the  following  goals: 

1 .  We  have  completed  the  construction  of  the  proposed  an  array  based  real-time  PAT  system, 
and  calibrated  and  tested  the  system  using  extensive  phantom  experiments. 

2.  We  have  built  the  animal  interface  and  successfully  tested  it  for  in  vivo  imaging  of  rat  brain. 

3.  We  have  conducted  phantom  experiments  that  confirmed  our  software  enhancement. 

Figs.  2a  and  2b  show  the  photograph  of  the  completed  PAT  system.  As  shown  in  Fig.  2a,  a 
tunable  Ti:Sapphire  laser  (part  A)  is  used  to  provide  690-950nm  near-  infrared  (NIR)  laser  pulses 
with  a  full-width  half-maximum  (FWHM)  value  of  8-25ns  pulse  width  with  a  repetition  rate  of  10Hz 
for  generating  photoacoustic  signals.  The  5  MHz  full-ring  transducer  array  (part  B)  is  used  for 
recording  the  photoacoustic  signals.  A  home  made  pre-  amplifier  (part  C)  receives  the  signals 
from  the  transducers  and  transmitted  the  amplified  signals  to  A/D  board  (part  D)  which  are 
then  stored  in  the  computer  (part  E). 
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Fig.  2.  Photograph  of  the  real-time  PAT  system,  (a):  Top  view,  (b):  Front  view. 


We  have  made  a  significant  progress  that  has  fulfilled  the  statement  of  work  proposed  for  Year  2 
of  this  project.  The  real-time  PAT  system  completed  in  Years  1  and  2  provides  us  a  platform  for 
performing  the  extensive  phantom  and  animal  experiments  proposed  for  the  coming  years. 

In  Task  3,  we  implemented  6  image  enhancement  schemes.  In  Year  1,  3  schemes  including  total 
variation  (TV),  advanced  regularization  techniques,  and  adaptive  meshing  were  developed 
successfully.  In  Year  2,  the  remaining  3  schemes,  including  dynamic  dual  meshing,  initial 
parameter  optimization,  and  reconstruction  of  absolute  optical  absorption  coefficient  for 
quantitative  high  resolution  functional  PAT  have  been  implemented  as  proposed. 


For  more  details  please  see  the  full  progress  report  of  the  second  year. 

3.  Software  Development  (Task  5) 

We  have  completed  the  software  development  proposed  for  Year  3  including  advanced 
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boundary  conditions,  3D  algorithms,  parallelization  of  3D  codes,  and  their  testing  using 
simulation  and  phantom  experiments. 


4.  Animal  Experiments  (Task  6) 

In  Year  3  we  conducted  extensive  animal  experiments  to  evaluate  the  PAT  system 
constructed  in  Years  1  and  2.  We  were  able  to  image  the  functional  anatomy  of  a  focal 
seizure  noninvasively  and  in  real  time  at  a  spatial  resolution  of  ~150pm  (Figs.  3  ). 

In  this  real  time  PAT  system  (Fig.  3a),  light  from  a  Ti:  Sapphire  laser  tunable  from  690  to  950  nm 
was  delivered  to  the  cortical  brain  surface  through  an  optical  fiber.  Part  of  the  energy  of  each 
laser  pulse  was  detected  by  a  photodiode  for  calibration.  A  192-element  full-ring  transducer 
array  was  used  to  capture  the  photoacoustic  signals  generated  by  the  laser  light.  All  the  192 
channels  were  equipped  with  preamplifiers  and  secondary  stage  amplifiers  for  optimized  signal- 
to-noise  ratio.  A  3:1  electronic  multiplexer  coupled  with  64-channel  analog-to-digital  converters 
completed  the  192-channel  data  acquisition,  providing  an  imaging  speed  of  0.33  s/frame  primarily 
limited  by  the  10  Hz  laser  repetition.  The  spatial  resolution  of  PAT  is  inversely  related  to  the 
bandwidth  of  the  ultrasonic  transducer.  In  our  PAT  system  each  ultrasonic  detector  had  a  5  MHz 
central  frequency,  a  70%  nominal  bandwidth,  and  a  diameter  of  6  mm  (Blatek,  Inc.,  PA,  USA), 
resulting  in  a  spatial  resolution  of  150pm. 

We  first  tested  the  PAT  system  accuracy  by  imaging  a  living  rat  brain  with  both  the  scalp  and  skull 
intact.  Results  indicate  that  the  noninvasive  PAT  image  of  the  rat  cortical  vasculature  correlated 
with  the  anatomical  location  of  each  vessel  (Figs.  38b,  c). 
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Figure  3  Real-time  PAT  system  for  seizure  dynamics,  (a)  Schematic  of  the  real-time  PAT  system.  A  192-element 
full-ring  transducer  array  was  used  to  capture  the  PA  signal  during  seizure  onset,  (b)  Noninvasive  PAT  imaging  of  a 
rat  brain  in  vivo  with  the  skin  and  skull  intact.  MCA,  middle  cerebral  artery;  RH,  right  hemispheres;  LH,  left 
hemispheres;  LOB,  left  olfactory  bulbs;  ROB,  Right  olfactory  bulbs,  (c)  Open-skull  photograph  of  the  rat  brain 
surface  acquired  after  the  PAT  experiment.  Numbers  1-5  indicate  the  corresponding  blood  vessels  in  the  PAT  image 
and  rat  brain  photograph. 

Accomplished  goals  in  the  third  year  of  the  project  are  summarized  as  follows: 

1.  We  have  developed  the  advanced  boundary  conditions  scheme,  and  implemented  the  3D 
reconstruction  codes  and  their  parallelization  which  allowed  efficient  high  performance  PAT 
image  reconstruction. 

2.  We  have  conducted  simulation  and  phantom  experiments  that  confirmed  our  software 
enhancement. 

3.  We  have  performed  extensive  in  vivo  experiments  using  the  rat  epilepsy  model  to  evaluate 
the  PAT  system  we  developed  in  this  project.  The  results  show  that  we  were  able  to  image 
the  functional  anatomy  of  epileptic  foci  noninvasively  and  in  real  time  at  a  spatial 
resolution  of  ~1 50pm. 

For  more  details  please  see  the  progress  report  of  the  third  year  of  the  project. 

5.  Rat  Model  of  Temporal  Lobe  Epilepsy  (Task  6) 

During  months  37-48  of  this  project,  we  have  successfully  developed  and  validated  the  rat 
model  of  temporal  lobe  epilepsy. 
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Accomplished  goals  in  the  fourth  year  of  the  project  are  summarized  as  follows: 


1 .  We  have  developed  and  validated  a  rat  model  of  temporal  lobe  epilepsy. 

2.  We  have  performed  extensive  in  vivo  experiments  to  achieve  concurrent  EEG  and  PAT 
recordings.  Interesting  results  are  obtained  that  may  be  significant  for  seizure  prediction  and 
for  neurovascular  coupling  studies. 

3.  We  have  designed  and  constructed  a  light-weighted  animal/transducer  interface  that  can  be 
used  for  chronic  monitoring  of  freely  moving  animals. 

After  many  trials  and  errors  in  Year  4,  we  have  successfully  fabricated  such  light-weighted 
transducer  array  (PVDF  film  transducers)  and  light  delivery  subsystem  (plastic  based  light 
delivery).  As  shown  in  Fig.  4,  a  cap  interface  firmly  attached  to  the  rat  head  will  allow  us  to 
chronically  monitor  seizure  activity  in  freely  moving  rat.  The  cap  interface  weights  only  ~7g  in 
total  (see  Fig.  10b).  In  Year  5  of  the  project  (i.e.,  12-month  no-cost-extension),  this  light- 
weighted  cap  interface  coupled  with  all  other  preparations  described  in  this  report  will  allow  us 
to  complete  the  proposed  Tasks  6  and  7. 
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Fig.  4  (a)  PAT/EEG  system/cap  interface  for  chronic  mo nitoring/im aging  of  freely  moving  rats,  (b) 

Schematic  detail  of  the  cap  interface 


For  more  details  please  see  the  progress  report  of  the  fourth  year  of  the  project. 


6.  Analysis  of  the  images  from  the  in  vivo  experiments  (Task  7) 

In  this  section,  the  accomplished  goals  in  the  year  5  of  the  project  are  presented  in  details. 

6.1.  Experimental  setup 

The  main  feature  of  the  photoacoustic  system  is  the  sparse  spherical  array  which  is  consisted  of 
192  discrete  transducers,  as  shown  in  Fig.  5.  The  transducers  were  mounted  on  a  custom 
fabricated  white  ABS  spherical  interface  with  a  distance  of  71 .6mm  from  the  center  of  the  array, 
and  formed  into  seven  layers.  Each  transducer  had  a  central  frequency  of  5MHz  and  a 
reception  bandwidth  of  greater  than  80%,  resulting  in  an  almost  isotropic  spatial  resolution 
about  0.2mm.  The  active  area  of  the  transducer  was  3mm  in  diameter,  and  the  angular 
acceptance  was  about  15  degree.  The  signal  from  the  192  transducers  was  amplified  by  16 
homemade  preamplifier  boards  and  coupled  into  a  64  channel  parallel  data  acquisition  system 
with  3:1  multiplexing,  as  indicated  with  three  different  colors. 
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Rats  were  elevated  to  the  center  of  the  spherical  interface  through  a  chamber  fixed  at  the  tank 
bottom,  whose  top  was  about  15  mm  beneath  the  interface  center,  and  a  transparent  plastic 
wrap  was  used  to  cover  the  chamber  top.  The  spherical  transducer  array  was  placed  in  the 
water  bath,  acoustically  coupled  to  the  rat  head  with  ultrasound  gel.  Laser  beam  was  delivered 
through  a  concave  lens  resulting  a  homogenously  illumination  on  the  rat  head.  In  this  study,  two 
wavelengths  of  1064nm  and  710nm  were  used.  The  1064nm  light  was  from  a  Q-switched 
Nd:YAG  laser,  and  the  710nm  light  was  from  a  tunable  Ti:sapphire  laser.  The  light  intensity  on 
the  rat  head  was  about  50mJ/cm2  for  1064nm,  and  4mJ/cm2  for  710nm,  all  below  the  ANSI 
exposure  limits  (ANSI  Z136).  Both  the  two  lasers  operate  at  10Hz,  so  that  one  complete  set  of 
3D  PAT  data  can  be  recorded  in  0.33s. 
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Fig.  5  Experimental  setup  of  the  combined  system.  The  anaesthetized  rat  was  elevated  to  the  center  of  the  spherical 
transducer  array  for  simultaneous  PAT  and  EEG  recording.  The  192-channel  signals  were  collected  and  multiplexed 
with  3:1  into  a  64  channel  DAQ  system,  and  with  a  10Hz  laser,  a  frame  rate  of  3.3f/s  is  achieved.  Then  the  collected 
data  is  stored  for  later  GPU  accelerated  reconstruction  using  3D  delay&sum  code  ; 

Two  syringe  needles  were  inserted  about  2mm  into  the  front  cortex  of  rat  brain  served  as  EEG 
electrodes.  One  electrode  was  chosen  as  reference  and  ground,  while  the  other  one  was  used 
as  signal  channel.  The  EEG  signal  was  amplified  (RA16PA,  TDtucker-Davis  Tech.)  and 
recorded  (RZ5  Bioamp  Processor,  TDtucker-Davis  Tech.)  with  a  sampling  rate  about  50kHz. 
Both  the  two  electrodes  were  firmly  glued  on  the  rat  head  and  the  EEG  cable  was  tightly  tied  on 
the  animal  holder  to  avoid  noise  induced  by  any  possible  unnecessary  movement  of  the 
operators.  The  PAT  system  gave  synchronize  signal  to  the  EEG  system  for  synchronization. 
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6.2  Animal  preparation 

Two  groups  of  rat  experiments  (n=16  rats,  about  35~45g,  with  8  rats  for  1064nm  and  the  other  8 
rats  for  710nm)  were  carried  out.  The  rats  were  firstly  anaesthetized  with  about  0.25ml  Urethane 
(0.25g/mL)  and  the  hair  on  the  top  head  was  shaven.  A  syringe  containing  0.04g/ml 
pentyienetetrazole  (PTZ)  with  a  scalp  needle  was  inserted  in  the  rat  abdomen  and  fixed  with 
glue.  This  was  to  easy  the  intraperitoneal  injection  later  in  the  experiment,  which  has  already 
been  established.  The  rat  was  immobilized  on  a  custom  made  plastic  holder  with  ear  bar  using 
glue  tapes.  Two  EEG  electrodes  were  inserted,  and  then  the  rat  was  elevated  into  the  center  of 
the  spherical  transducer  array.  The  seizure  process  was  recorded  for  25min,  and  PTZ  was 
administrated  with  a  dose  of  0.5mL  at  5min.  The  injection  of  the  drug  lasted  about  30s.  All  rats 
kept  alive  through  the  whole  data  recording  process.  All  the  animal  procedures  performed  were 
in  accordance  with  the  approved  University  of  Florida  lACUC  protocols. 

6.3  Data  processing 

With  a  frame  rate  of  3.3f/s,  there  were  5000  complete  sets  of  PAT  data  for  the  25min 
experimental  time.  Before  doing  3D  reconstruction  and  processing,  2D  time  series  images  were 
reconstructed  with  the  signals  from  the  64  transducers  in  the  lowest  layer  of  the  array,  which 
were  indicated  with  blue  in  Fig. 5.  Then  all  2D  time  series  were  first  screened  for  movement 
artifacts  using  a  movie  function,  although  rats  were  completely  paralyzed  during  experiments.  In 
our  experiments,  only  3  out  of  8  rats  for  the  1064nm  didn’t  satisfy  this  condition,  while  all  the  8 
rats  for  710nm  were  quantified.  Then  3D  images  were  reconstructed  using  3D  delay&sum 
algorithm,  which  was  accelerated  by  GPU  parallel  technique.  The  resulted  3D  images  had  a 
voxel  number  of  201x201x101  representing  a  10x10x5mm  volume,  and  the  negative  voxel 
values  were  set  to  zero.  The  reconstruction  speed  is  about  0.97s  for  one  frame.  Amira  (from 
Visage  Imaging,  Inc.)  was  used  to  render  the  3D  images. 

There  were  two  sources  of  noise  in  the  EEG  signal,  a  strong  sharp  10Hz  noise  from  the  Q- 
switch  of  the  laser,  and  the  60Hz  alternating  current  noise.  For  de-noising,  the  data  deteriorated 
by  the  Q-switch  was  replaced  with  adjacent  data,  and  the  60Hz  alternating  current  noise  is 
removed  by  filtering  out  58-62Hz  in  frequency  domain  (Fig.  6  b).  Then  the  EEG  data  was 
applied  with  a  low  pass  filter  of  100Hz,  and  transformed  into  a  300000  point  signal  for  the  25min 
EEG  recording,  corresponding  to  sampling  rate  of  200Hz. 

6.4  Neurovascular  coupling  study 

Because  the  1064nm  and  710nm  light  were  mainly  absorbed  by  oxy-hemoglobin  (Hb02)  and 
deoxy-hemoglobin  (Hbr)  in  the  tissue  respectively,  with  reasonable  approximations  the  PAT 
signals  from  them  can  be  taken  as  the  total  volume  change  of  those  two  absorbers.  For 
neurovascular  coupling  study,  averaged  signals  from  specific  anatomical  regions  of  interest 
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(ROI)  are  compared  to  that  from  the  superior  sagittal  sinus  (SSS)  as  well  as  the  EEG  signal. 
The  PAT  signals  were  expressed  as  the  percentage  changes  from  the  baseline,  which  is 
defined  to  be  the  average  value  over  first  4min  of  the  recording.  Next,  we  further  investigated 
the  correlation  relationship  between  these  signals.  We  also  calculated  the  correlation  coefficient 
of  every  voxel  with  the  SSS,  and  got  the  negative  correlation  map  following  the  seizure  onset  to 
identify  important  regions  in  the  brain.  We  choose  the  PAT  signal  from  the  SSS  because  it’s  well 
distinguished  in  the  reconstructed  3D  images  and  has  the  highest  SNR.  In  this  step,  the 
reconstructed  3D  image  was  applied  with  a  moving  window  of  9-voxel  wide  cube  for  averaging 
to  improve  the  low  SNR  due  to  the  limited  number  of  transducers.  Finally,  the  granger  causality 
networks  between  the  SSS,  hippocampus,  and  cortex  at  different  times  were  calculated.  All  the 
acquired  results  were  examined  by  experts  in  the  neurology. 

6.5  Results 

The  EEG  signal  was  stable  before  the  drug  injection.  Ictal  discharges  were  induced  shortly  after 
the  injection  of  PTZ,  which  were  characterized  by  the  low-amplitude  fast  activity  with  high 
amplitude  (Fig.  6  A).  The  epileptic  pattern  progressively  increased  in  amplitude  in  the  following  2 
to  3  minutes.  As  the  seizure  evolves,  it  gradually  turned  into  clonic  paroxysmal  discharges  of 
lower  frequency.  2D  plot  of  the  time-frequency  power  map  of  the  EEG  signal  shows  that  a 
mixture  of  lower  signal  frequency  particularly  in  1-3  Hz  was  seen  throughout  the  seizure,  and 


higher  frequency  power  components  (up  to  ~4-6  Hz)  were  stronger  early  in  the  seizure  and 


gradually  decreased  in  power  (Fig.  6  B).  This  is  close  to  results  in  other  publications.  Seizures 
were  induced  successfully  in  all  the  16  rats,  and  they  were  similar  in  the  shape  and  amplitude. 
The  injection  time  is  easily  distinguished  as  isolated  high  amplitude  spikes  around  5min  in  the 
EEG  signal. 
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Fig.  6  Electrophysiology  of  PTZ  induced  seizure,  (a)  An  example  of  EEG  traces.  The  arrow  highlights  the  time  of  PTZ 
injection.  The  seizure  always  starts  with  synchronized  high-frequency  oscillations  that  increase  in  amplitude, 
progressing  into  clonic  paroxysmal  discharges  of  lower  frequency  afterwards,  (b)  Time-frequency  analysis  for  the 
same  EEG  signal  above.  Analysis  was  performed  using  a  short-time  Fourier  transform.  Peak  power  was  normalized 
to  1.  The  dominant  ictal  signal  frequency  was  1-3  Hz.  Higher  frequency  components  up  to  6  Hz  occurred  earlier  in 
seizure  and  diminished  in  power  gradually. 

6.6  PAT  image  reconstruction 

Representative  3D  images  for  the  Hb02  and  Hbr  are  displayed  in  Fig.  7A  and  B  compared  with 
a  photo  of  the  rat  brain  in  Fig.  7  C.  Although  the  image  quality  is  not  excellent  due  to  the  limited 
transducer  numbers,  the  SSS  and  the  transverse  sinuses  (TS)  are  clearly  reconstructed.  It’s 
also  noted  that  in  the  3D  image  of  Hb02  the  Ostia  end  of  the  straight  sinus  (SS)  which  is  below 
the  SSS  in  z  direction  can  also  be  seen.  The  SNR  of  the  3D  image  for  Hbr  is  not  as  good  as  that 
for  Hb02  partly  due  to  relatively  low  laser  intensity,  and  the  other  reason  lies  in  the  fact  that  the 
concentration  of  Hb  is  much  lower  than  that  of  Hb02  in  the  blood  vessels. 
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Fig.  7  The  reconstruction  of  3D  PAT  images.  (A)  and  (B)  show  the  representative  3D  images  for  Hb02  and  Hbr  in 
different  views  respectively.  The  image  domain  is  10x10x5mm  containing  201x201x101  pixels.  Thresholds  and 
colorbars  are  presented  for  the  3D  images,  (b)  Photography  of  rat  head  with  hair  removed.  Scale  bar  represents  2mm. 
The  SSS  and  TS  are  well  distinguished  in  both  the  Hb02  and  Hbr  images,  and  the  Ostia  end  of  the  inferior  sagittal 
sinus  can  also  be  seen  in  the  Hb02  image. 

6.7  Neurovascular  coupling  in  selective  regions 

To  investigate  the  hymodynamic  changes  during  the  seizure  initiation  and  evolution,  regions  of 
interest  (ROI)  are  selected  for  the  SSS,  as  well  as  other  regions  near  the  hippocampus  and  the 
cortex  in  the  brain  that  are  reported  to  have  high  gamma-aminobutyric  acid  (GABA)  receptor 
concentrations  (Fig.  8  and  9).  ROI  was  also  chosen  in  the  upper  posterior  region  in  the 
reconstruction  domain,  which  contained  no  tissue  of  the  rat  served  as  reference.  With  a  high 
frame  rate  of  3.3f/s,  the  rapid  hemodynamic  changes  in  these  regions  \can  be  captured 
photoacoustically.  Taking  the  SSS  and  TS  as  positioning  references,  it’s  easy  to  overlay  the 
maximum  amplitude  projection  (MAP)  PAT  vascular  structural  images  precisely  onto  the  high 
resolution  anatomical  MRI  slices  of  the  rat  brain  obtained  from  online  database  . 

Fig.  8  shows  the  averaged  Hb02  signals  of  the  ROIs  during  the  seizure,  with  corresponding 
EEG  signal  presented  to  confirm  the  onset,  duration,  and  morphology  of  the  seizure.  The  Hb02 
signals  are  converted  into  the  relative  percentage  change  to  the  baseline.  Because  the  SSS  has 
the  highest  SNR  in  the  reconstructed  3D  images,  the  SD  for  SSS  is  much  higher  than  the  other 
three  signals.  The  Hb02  signals  and  the  EEG  signal  are  mutually  well  correlated.  All  the  Hb02 
signals  of  the  SSS,  hippocampus  and  cortex  remained  quite  still  at  first.  We  can  largely  divide 
the  change  of  these  three  signals  after  the  injection  of  PTZ  into  three  for  stages,  which  is 
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indicated  in  Fig.  8  B.  We  can  see  that  there  was  an  ‘epileptic  dip’  in  the  SSS  at  the  seizure 
initiation,  and  a  prominent  increase  of  Hb02  was  shown  in  the  hippocampus  at  the  same  time 
which  was  probably  induced  by  the  dilation  of  the  arteries.  It’s  also  seen  that  the  change  of 
hippocampus  is  largely  negatively  correlated  with  the  SSS,  while  the  cortex  and  the  SSS  are 
positive  correlated.  In  the  third  stage,  both  the  Hb02  and  Hbr  signals  in  the  SSS  and  cortex 
dropped,  while  the  increase  of  Hb02  in  the  hippocampus  sustained  through  the  whole 
experiment.  Comparably,  the  reference  signal  remains  the  same  all  although  the  experiment. 


baseline  1st  2nd  3rd 


Fig.  8  The  Hb02  changes  of  the  SSS,  hippocampus  and  cortex.  (A)  3D  vascular  images  of  the  rat  brain  (a1)  and  the 
MAP  vascular  structural  images  (a2  and  a3  in  the  x-y  and  x-z  plane  respectively)  overlaid  on  the  MRI  anatomical 
images  using  1064nm.  The  ROIs  for  the  SSS,  hippocampus,  cortex,  and  the  reference  are  indicated.  Thresholds, 
colorbars,  and  scale  bar  are  presented.  (B)  The  Hb02  changes  of  the  SSS,  hippocampus,  cortex  and  the  reference 
as  well  as  the  corresponding  EEG  signal.  The  changes  of  the  Hb02  signals  are  expressed  in  percentage  relative  to 
the  baseline,  which  is  defined  as  the  first  four  minutes  of  the  recording.  The  changes  of  these  Hb02  signals  after  the 
PTZ  injection  are  divided  into  three  stages,  as  indicated. 
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Fig.  9  presents  result  in  another  rat  showing  how  the  changes  of  Hbr  in  the  SSS,  hippocampus, 
cortex  and  the  EEG  signal  are  mutually  correlated.  The  sizes  and  the  relative  positions  of  the 
three  ROIs  are  defined  to  be  the  same  as  those  in  Fig.  8A.  The  EEG  signal  here  looked  almost 
the  same  in  pattern  with  that  in  Fig.  8  B,  indicating  that  our  PTZ  induced  generalized  seizure 
modal  was  quite  stable  and  the  drug  administration  was  well  controlled,  so  that  results  in  Fig.  8 
and  Fig.  9  were  comparable.  We  also  divided  the  timecourse  of  the  Hbr  changes  into  the  three 
stages  with  the  exact  timing  as  that  in  Fig.  8  for  easy  comparison.  Although  with  a  relatively  low 
SNR,  we  clearly  see  that  the  decrease  of  Hbr  in  the  SSS  lagged  the  notable  ‘epileptic  dip’  in  the 
SSS  shown  in  Fig.  8  B.  In  addition,  there  was  a  slightly  decrease  of  Hbr  in  the  hippocampus 
followed  by  an  increase  in  the  first  and  second  stages,  but  these  changes  were  not  as  notable 
as  the  increase  of  Hb02  in  the  hippocampus  at  the  same  time  in  Fig.  8  B.  In  general,  the  Hbr 
signals  for  the  hippocampus  and  cortex  were  generally  the  same  in  shape,  which  were 
negatively  correlated  with  the  Hbr  of  SSS. 


Fig.  9  The  Hbr  changes  of  the  SSS,  hippocampus  and  cortex.  (A)  3D  vascuiar  images  of  the  rat  brain  (a1)  and  the 
MAP  vascular  structural  images  (a2  and  a3  in  the  x-y  and  x-z  plane  respectively)  overlaid  on  the  MRI  anatomical 
images  using  710nm.  The  ROIs  for  the  SSS,  hippocampus,  cortex,  and  the  reference  are  indicated.  Thresholds, 
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colorbars,  and  scale  bar  are  presented.  (B)  The  Hb02  changes  of  the  SSS,  hippocampus,  cortex  and  the  reference 
as  well  as  the  corresponding  EEG  signal.  The  changes  of  the  Hbr  signals  are  expressed  in  percentage  relative  to  the 
baseline,  which  is  defined  as  the  first  four  minutes  of  the  recording.  The  changes  of  these  Hb02  signals  after  the  PTZ 
injection  are  also  divided  into  three  stages  with  the  same  timing  as  that  in  Fig.  4B,  as  indicated. 

Fig.  10  shows  the  statistical  results  of  the  Hb02  and  Hbr  changes  across  all  the  rat  experiments, 
which  are  similar  to  those  in  Fig.  8  and  9.  The  average  signal  in  each  rat  is  calculated  every  30 
seconds,  and  then  further  averaged  between  different  rats,  and  the  corresponding  SD  is 
calculated  and  shown. 


Fig.  10  The  average  Hb02  and  Hbr  changes  of  the  SSS,  hippocampus  and  cortex  across  all  the  rats.  The  signals  are 
averaged  every  30  seconds,  and  corresponding  SD  between  different  rats  are  also  calculated  and  shown  (in  red).  (A) 
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and  (B),  Hb02  and  Hbr  changes  of  the  SSS;  (C)  and  (D)  Hb02  and  Hbr  changes  of  the  hippocampus;  (E)  and  (F) 
Hb02  and  Hbr  changes  of  the  cortex. 

6.8  Correlation  analysis  and  hyper-oxygenation  localization 

Based  on  the  above  results,  the  hippocampus  showed  strong  hyperfusion  and  hyperoxgenation 
during  the  seizure,  and  the  Hb02  change  of  it  was  inversely  related  with  that  of  the  SSS.  Fig.  1 1 
A  shows  the  time  change  of  the  correlation  coefficient  between  the  Hb02  signals  of  these  two 
regions,  and  it’s  found  that  this  value  was  positive  before  the  drug  injection  (Fig.  1 1 B),  but  it 
turned  negative  for  most  of  the  time  after  the  seizure  onset  (Fig.  1 1 C).  All  the  correlation 
coefficients  in  this  paper  were  calculated  in  a  short  time  window  of  1  min. 


Fig.  11  Correlation  study  of  the  Hb02  signals.  (A)  The  time  change  of  the  correlation  coefficient  between  the 
hippocampus  and  SSS.  (B)  and  (C)  the  Hb02  signal  changes  of  the  SSS,  hippocampus  and  cortex  at  two  different 
time  window  as  indicated  in  A  with  blue  rectangles.  (D)  The  3D  negative  correlation  coefficient  map  at  the  time  of  the 
green  vertical  bar  in  A,  shown  in  different  views,  and  the  MAP  images  of  it  in  the  x-y  and  x-z  plane  imposed  on  the 
corresponding  MAP  Hb02  vascular  images  and  the  MRI  anatomical  images.  (E)  The  3D  negative  correlation 
coefficient  map  at  different  times  as  indicated  in  A  with  vertical  bars.  The  threshold  for  the  negative  correlation  map  is 
chosen  to  be  0.6  in  D,  and  0.4  in  E  and  F.  It’s  seen  that  the  negative  correlation  voxel  with  the  SSS  for  Hb02  didn’t 
shown  before  the  drug  injection,  and  were  concentrated  in  the  hippocampus  and  its  surroundings  after  the  seizure 
initiation. 
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For  further  investigation,  we  calculated  the  3D  Hb02  negative  correlation  coefficient  map 
between  every  voxel  and  the  SSS  (Fig.  1 1  D  and  E).  In  this  step,  the  3D  Hb02  images  were 
firstly  smoothed  with  a  moving  window  of  a  9-voxel  wide  cube  to  improve  the  SNR.  Fig.  11  D 
presents  the  3D  negative  correlation  coefficient  map  at  different  times,  as  indicated  in  Fig.  1 1  A 
with  vertical  red  lines.  Fig.  11  E  shows  the  same  result  as  Fig.  7  D  (e),  but  from  different  views, 
along  with  the  MAP  images  of  it  in  the  x-y  and  x-z  planes  coregistered  with  the  corresponding 
MAP  Hb02  vascular  structural  images,  and  the  MRI  anatomical  images.  We  found  that  voxel 
with  high  negative  coefficient  (with  a  threshold  of  0.7)  were  hardly  found  before  the  drug 
injection,  but  they  began  to  show  and  concentrated  in  regions  around  the  hippocampus  after  the 
seizure  onset,  although  asymmetry  is  present  in  these  images  (Fig.  11E).  Especially,  when  the 
seizure  just  began  later  in  the  first  stage  of  the  seizure  (around  6-7min),  the  Hb02  signal  in  the 
hippocampus  was  correlated  with  that  of  the  SSS  with  a  large  negative  correlation  coefficient, 
so  that  the  3D  negative  correlation  map  is  preventative  at  this  time  (Fig.  11  D).  We  also 
calculated  the  negative  correlation  maps  of  Hb02  in  other  rat  experiments,  and  similar  results 
were  obtained  (Fig.  12).  In  summery,  voxel  with  high  negative  correlation  coefficient  began  to 
show  after  the  seizure  initiation  in  the  hippocampus  and  its  surroundings,  which  indicates  that 
the  hippocampus  seemed  to  be  more  sensitive  in  this  intraperitoneal  PTZ  seizure  model.  This 
also  implies  that  this  method  may  be  of  great  help  for  seizure  localization  especially  in 
generalized  seizure  cases  where  fMRI  is  of  less  help. 


Fig.  12  Negative  correlation  map  of  Hb02  in  different  rats.  The  time  window  for  generating  the  correiation  map  is 
chosen  to  be  from  6min  to  7min.  The  results  are  presented  in  2D  maps,  which  are  imposed  on  the  corresponding 
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MAP  Hb02  vascular  images  and  the  MRI  anatomical  images  in  x-y  and  x-z  planes,  and  also  presented  in  3D  images, 
which  are  coregistered  with  the  3D  vascular  images.  The  threshold  for  the  negative  correlation  map  is  chosen  to  be 
0.6  in  the  2D  images,  and  0.4  in  the  3D  images.  (A)-(D)  represent  the  results  for  different  rats. 

6.9  Granger  causality  analysis 

In  Fig.  11  B  and  C,  we  see  that  the  hemodynamic  coupling  of  the  high  frequency  fluctuations 
between  different  brain  regions  is  different  with  different  times.  This  implies  that  the  high  ‘noise’ 
component  of  the  PAT  signals  actually  contains  rich  seizure  relevant  information  that  can  not  be 
inferred  from  the  overall  trend.  We  applied  short  time  window  spectral  pairwise  granger 
causality  method  to  identify  the  functional  interaction  among  the  SSS,  hippocampus  and  cortex 
during  the  seizure  process.  Granger  causality  allowed  for  assessment  of  the  magnitude  and 
direction  of  temporal  relationships  among  multiple  signals  during  overlapping  time-series 
windows,  holding  the  promise  of  revealing  paths  of  information  flow  within  the  nervous  system. 
This  method  has  emerged  as  the  leading  statistical  quantities  to  furnish  directional  information 
from  multivariate  neural  data. 

To  make  use  of  this  formulation,  the  PAT  signals  in  each  short  time  window  is  firstly  detrend 
with  a  3'^'^  order  and  then  normalized  by  dividing  the  standard  deviation.  In  this  way,  the  PAT 
signals  can  be  approximately  regarded  as  a  stationary  stochastic  process  with  locally  stationary 
segments,  and  the  collection  of  the  PAT  signals  from  different  rat  experiments  can  be  treated  as 
an  ensemble  of  realizations.  In  this  paper,  the  time  window  is  set  to  be  Imin,  which  is  200  points 
for  the  frame  rate  of  3.33Hz.  This  is  a  promise  between  preserving  the  time  variability  and 
maintaining  the  smoothness  of  the  estimated  spectral  quantities.  The  model  order  is  determined 
with  theAkaike  Information  Criterion  (AlC).  Because  the  AlC  generally  decreases  monotonically 
with  increasing  order  for  all  the  windows,  we  set  the  model  order  to  be  19  selected  as  a  tradeoff 
between  sufficient  spectral  resolution  and  overparameterization,  and  the  model  order  is 
validated  with  white  noise  test.  The  directional  causality  was  calculated  as  the  average  value  in 
the  0-0. 3Hz  range  which  makes  up  about  84%  of  the  total  power  in  the  1.6Hz  frequency  band. 
Granger  causality  spectral  analysis  of  both  the  Hb02  and  Hbr  data  were  performed  for  all  pairs 
between  the  SSS,  hippocampus,  cortex,  and  the  reference  signals  which  is  used  for  significance 
assessments. 
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Fig.  13  Granger  causality  analysis  of  the  Hb02  and  Hbr  signals  between  different  regions  of  the  brain  and  the 
reference  signals.  (A)  and  (B)  the  granger  causality  of  Hb02  between  the  SSS,  hippocampus  and  cortex;  (C)  and  (D) 
the  granger  causality  of  Hb02  between  the  above  three  brain  regions  and  the  reference;  (E)  and  (F)  the  granger 
causality  of  Hbr  between  the  SSS,  hippocampus  and  cortex.  The  horizontal  dotted  lines  in  (A)  and  (B)  are  the 
average  value  of  the  granger  causality  of  Hb02  between  these  three  brain  regions  and  the  reference,  or  rather  the 
average  value  of  (C)  and  (D).  Similarly,  the  horizontal  dotted  lines  in  (E)  and  (F)  are  the  average  value  of  the  granger 
causality  of  Hbr  between  these  three  brain  regions  with  the  reference.  SSS,  superior  sagittal  sinus;  Hipp, 
hippocampus;  Cort,  cortex. 
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Fig.  13  A  and  B  show  the  time  varied  granger  causality  of  Hb02  between  pairs  among  the  SSS, 
hippocampus  and  cortex.  Fig.  13  C  and  D  are  the  granger  causality  of  Hb02  between  these 
three  brain  regions  with  the  reference,  which  is  less  varied  with  time,  and  the  average  value  of  it 
can  be  taken  as  a  noise  background  or  baseline,  as  indicated  in  Fig.  13  A  and  B  with  dotted 
black  horizontal  lines.  Similarly  Fig.  13  E  and  F  show  the  granger  causality  of  Hbr  for  these 
three  brain  regions.  The  results  of  Hb02  shows  that  strong  influence  of  the  hippocampus  to  the 
SSS  is  present  in  the  first  and  second  stages  of  the  seizure;  strong  influence  of  the  cortex  to  the 
SSS  is  seen  in  the  second  and  early  third  stages;  and  as  the  seizure  evolves,  significant 
causality  for  all  pairs  between  these  three  regions  can  be  seen  in  both  directions  late  in  the  third 
stage.  Comparatively,  the  causality  of  the  Hbr  for  pairs  between  these  three  brain  regions  is 
overall  less  significant,  which  may  be  due  to  the  poor  SNR  of  the  data.  However,  some  time 
relevant  features  can  still  be  obtained  through  the  calculated  results. 

The  granger  causality  in  different  times  can  be  better  presented  in  graphs,  as  shown  in  Fig.  14. 
In  Fig.  14,  The  arrow  presents  the  direction  of  the  influence,  and  the  visibility  presents  the 
magnitude  of  the  causality.  The  granger  causality  of  less  significance  is  not  shown,  and  the 
thresholds  for  Hb02  and  Hbr  are  chosen  to  be  the  average  value  of  the  causality  between  the 
reference  and  the  reference,  or  rather  the  dotted  lines  in  Fig.  13  A  (or  B)  and  E  (or  F) 
respectively.  The  value  of  the  causality  can  be  read  from  the  colorbars.  We  see  that  as  the 
seizure  evolved,  the  connectivity  between  these  three  brain  regions  tended  to  be  more  and 
much  closer. 


A 
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Fig.  14  Granger  causality  graphs  of  Hb02  (A)  and  Hbr  (B)  between  different  regions  of  the  brain.  The  granger 
causality  values  are  coded  by  the  visibility  of  the  lines  between  these  brain  regions  (Center),  and  the  arrowheads 
indicate  the  direction  of  granger  causal  influence.  Lines  between  site  pairs  not  reaching  significance  in  the  coherence 
and  Granger  causality  measures  are  not  shown.  SSS,  superior  sagittal  sinus;  Hipp,  hippocampus;  Cort,  cortex. 

6.10  Significance 

Epilepsy  affects  about  1%  of  the  population  world  wide,  and  about  20  to  30%  of  these  patients 
are  refractory  to  all  forms  of  medical  treatment.  For  those  patients  the  best  treatment  option  is 
resective  brain  surgery,  which  however  critically  depends  on  the  complete  mapping  of  the  local 
epileptic  seizures,  and  the  well  understanding  of  the  possible  epileptic  circuit. 

Much  of  our  understanding  about  seizure  comes  from  electrophysiological  recording  methods, 
and  intracranial  EEG  is  commonly  used  for  presurgery  seizure  localization.  However,  its  defects 
of  invasiveness,  low  spatial  sampling  and  volume  conduction  may  leads  to  incomplete  seizure 
mapping.  As  a  result,  various  neuroimaging  technologies  have  been  employed  to  give  a  high- 
resolution  real  time  spatial  and  temporal  ‘read  out’  of  the  dynamics  of  cortical  processing. 

As  the  most  common  neuroimaging  method,  fMRI  is  well  applied  for  the  diagnosis  and 
presurgery  localization  of  focused  seizure,  where  the  hemodynamic  changes  mainly  occurred 
within  the  seizure  focus  and  can  be  well  distinguished  in  fMRI.  However,  fMRI  can  not  localize 
generalized  seizure  because  in  this  case  the  hemodynamic  changes  can  be  found  wide-spread 
in  the  brain,  and  it  remains  unclear  whether  they  arise  from  a  portion  of  the  brain  or  from  more 
widespread  circuits  . 

This  predicament  mainly  owes  to  two  reasons.  One  reason  is  that  the  spatial  and  temporal 
resolution  of  fMRI  is  poor,  so  the  initial  site  of  seizure  is  missed.  Most  importantly,  fMRI  is  only 
sensitive  to  Hbr.  It  relies  on  the  drop  of  Hbr  in  the  vein  due  to  the  speedy  of  the  blood  flow  to 
localize  the  activation  sites  of  the  neurons.  However,  the  change  of  Hbr  and  the  blood  flow  is 
two  counter  parts  during  the  seizure,  so  that  the  changes  of  Hbr  in  the  tissue  may  not  be  as 
notable  as  that  of  Hb02,  which  mainly  results  from  the  dilation  of  the  artery  and  the  blood  flow. 

In  generalized  seizure,  although  current  diagnostic  methods  including  intracranial  EEG  and 
fMRI  all  find  that  regions  all  over  the  brain  are  involved  in  seizure,  this  doesn’t  surely  means  that 
all  these  regions  are  the  origin  of  seizure  or  equally  in  function.  Primary  regions  or  seizure 
network  may  be  revealed  if  more  information  is  given.  With  the  realtime  PAT  operated  at  two 
different  wavelengths,  relative  focused  regions  that  were  negatively  correlated  with  the  SSS 
were  found  in  the  hippocampus  and  its  surroundings  following  the  onset  of  the  seizure  using  the 
Hb02  data,  and  the  change  of  Hb02  in  the  hippocampus  was  found  to  be  ahead  of  the  Hbr  in 
the  PTZ  induced  generalized  seizure  model  (Fig.  8,  9,  and  10).  Since  the  SSS  is  decreased  in 
Hb02  when  the  seizure  began,  the  above  regions  with  high  negative  correlation  coefficient 
represent  regions  with  hyperoxygenation  and  hyperfusion,  which  is  usually  the  marker  of  the 
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seizure  location.  PTZ  is  a  GABA  antagonist,  and  this  result  is  in  coincidence  with  the  findings 
hippocampus  is  high  in  GABA  receptor  concentration.  All  these  imply  that  PAT  may  serve  better 
for  seizure  localization.with  appropriate  laser  wavelength  settings  and  data  processing  methods. 

Besides,  generalized  seizure  occurs  in  multiple  regions  all  over  the  brain  that  form  a  seizure 
network.  It’s  important  to  know  the  relationship  and  function  of  these  involved  regions.  The 
spatial  and  temporal  resolutions  of  our  realtime  PAT  system  are  both  about  10  times  higher  than 
that  of  fMRI,  thus  more  subtle  hymodynamic  changes  related  to  the  seizure  could  be  identified. 
We  applied  short  time  window  spectral  pairwise  granger  causality  method  to  both  the  Hb02  and 
Hbr  data,  and  the  time  resolved  interaction  between  different  regions  of  the  brain  was  visualized, 
resulting  in  a  better  understanding  of  epileptic  networks  at  high  spatiotemporal  resolution. 

6.11  Data  processing  method 

PAT  has  already  been  employed  in  the  neurovascular  coupling  study  with  various  animal 
models,  which  however  are  mainly  limited  in  the  monitoring  of  hemodynamic  signals  in  large 
blood  vessels.  Actually  the  changes  in  capillaries  are  more  related  to  the  metabolic  status  of  the 
cells.  Yet  the  direct  monitoring  of  the  capillary  is  difficult  because  the  resolution  of  PAT  under 
the  epidermis  is  usually  not  enough  due  to  the  diffusion  of  the  light.  In  this  work,  the  signals  are 
averaged  in  the  surrounding  region  to  increase  the  SNR,  and  converted  into  percentage 
changes  regarding  to  the  baseline  for  comparison.  With  this  method,  the  small  changes  in  deep 
tissue  that  is  submerged  in  the  signal  of  larger  blood  vessels  can  be  revealed.  For  example.  Fig. 
15  shows  a  2  minute  signal  summed  within  a  9  voxel  wide  cube  around  the  hippocampus  with 
1064nm  (the  black  line  in  the  upper  right  corner),  expressed  as 

h  =  '^Aj,k  =  X  ^ij,k  +  X  ^J,k  =  A  +  ^2 

mod(/+ j+k,2)=0  mod(/+ j+k,2)=l 

compared  with  /j  and  1^  (lines  in  the  lower  right  corner  in  red  and  blue  respectively,  which  are 
the  summed  signals  with  the  two  halves  of  the  voxel  in  the  cube),  where  value  of  a 

voxel  in  the  cube.  The  difference  of  /j  and  I^can  almost  be  neglected,  and  calculation  shows 
that  the  SD  of  /j  -1^  is  just  0.04  of  the  SD  of  1^.  This  means  the  white  noise  can  almost  be 
neglected.  Actually,  the  high  frequency  fluctuations  in  the  PAT  data  here  are  more  related  to  the 
hemodynamic  change  in  the  tissue,  which  can  be  inferred  from  the  correlation  between  the 
hippocampus  and  the  SSS.  The  Hb02  signals  in  the  hippocampus  and  SSS  are  positively 
correlated  before  seizure,  but  negatively  correlated  with  a  high  correlation  coefficient  shortly 
after  the  seizure  onset.  This  can  be  totally  explained  by  the  fluctuation  of  the  laser  or  the 
influence  of  the  reconstruction  algorithm. 


22 


Huabei  Jiang,  PhD/Final  Report 


.  7  77  7i  36  76  3  33  34  36  38  4 

dmestimn) 


Fig.  15  White  noise  assessment  of  the  reconstructed  images.  Upper  right  corner  shows  the  summed  signal  of  a  9 
voxel  wide  cube,  compared  with  the  signals  summed  with  the  two  halves  of  the  total  voxel  in  the  cube  (the  red  and 
blue  lines  in  lower  right  corner). 

6.12  Other  factors  in  neurovascular  coupling 

Contradictory  data  exists  for  neurovascular  coupling.  One  reason  for  this  lies  in  the  model 
difference  and  the  individual  difference,  as  well  as  the  systematic  influence  of  the  animals  ,  and 
the  other  reason  may  caused  by  the  different  principles  of  various  technologies  leading  to 
different  conclusion  in  case  of  incomplete  acquisition  of  the  information. 

Neurovascular  coupling  refers  to  all  the  aspects  regarding  the  neural  metabolic  demands  and 
the  energy  source  supply  by  the  vascular  system,  including  the  local  filed  potential,  the  oxygen- 
consumption  rate,  oxygen  partial  pressure,  Hb02  and  Hbr  volume,  blood  flow,  and  so  on.  More 
parameters  we  can  measure,  the  more  correct  understanding  of  the  neurovascular  coupling  we 
can  get  for  the  diagnosis  of  epilepsy. 

In  this  work,  through  the  measuring  of  Hb02  and  Hbr,  other  factors  such  as  the  dilation  of  the 
arteries  and  blood  flow,  which  works  as  means  of  power  supply  to  the  neurons  from  the 
vascular  system,  can  also  be  indirectly  inferred.  However,  signal  analysis  can  still  be 
complicated  sometimes.  We  only  get  information  of  Hb02  and  Hbr,  making  it  difficult  to  infer  the 
status  of  the  neurons.  For  example,  we  observed  that  the  Hb02  is  decreasing  while  the  Hbr  is 
increasing  at  the  second  stage  in  hippocampus,  but  since  the  overall  blood  supply  is  increase  in 
the  rat  brain,  we  are  not  sure  whether  the  temporal  decreasing  of  the  Hb02  in  the  hippocampus 
stands  for  the  decreasing  of  neural  activity  there  or  it’s  caused  by  insufficient  blood  supply,  and 
which  factor  dominates.  Furthermore,  the  oxygen  partial  pressure  in  the  tissue  which  is  more 
related  to  the  neural  activity  can  be  uncoupled  with  the  Hb02  change.  This  also  impedes  our 
understanding  of  neurovascular  coupling. 

The  preceding  of  the  hemodynamic  signal  to  the  EEG  signal  in  ictal  seizures  has  also  been 
reported  with  various  kinds  of  technologies.  However,  in  other  cases  and  in  our  studies,  this 
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phenomenon  is  not  seen.  Because  this  is  of  quite  importance  in  aspect  of  seizure  prediction  and 
manual  intervention,  more  ictal  neurovascular  coupling  studies  should  be  done  with  different 
seizure  models  to  find  out  the  condition  and  the  mechanism  of  this  phenomenon. 

6.13  Limitation  in  our  work 

There  are  limits  in  our  work.  First  of  all,  the  image  contrast  and  resolution  are  greatly  hindered 
by  the  limited  transducer  numbers  with  our  system.  Besides,  although  the  morphology  of 
seizures  exhibit  quite  similarity  for  all  the  rats  in  our  experiments,  images  with  the  two  different 
wavelengths  are  not  acquired  on  the  same  animal.  Although  this  will  be  not  likely  to  refute  the 
main  conclusions  we  got,  valuable  information  that  need  simultaneous  dual-wavelength 
recording  to  be  unearthed  is  missed.  Finally  and  most  important  of  all,  generalized  seizure 
involves  wide-spread  network  all  over  the  brain,  in  which  the  electrical  and  hemodynamic  signal 
of  different  regions  are  closely  related.  In  our  work,  we  only  analyzed  the  granger  causality  of 
the  hemodynamic  change  in  SSS,  hippocampus,  and  cortex.  Detailed  study  needs  to  be  done  to 
dig  out  the  interaction  of  more  different  regions. 

7.  Simultaneous  real-time  3D  photoacoustic  tomography  and  EEG  for  neurovascular 
coupling  study  in  an  animal  model  of  epilepsy 

A.  System  setup 

The  real-time  3D  photoacoustic  system  has  already  been  reported,  and  the  main  feature  of  this 
system  is  the  sparse  spherical  array  which  is  consisted  of  192  discrete  transducers,  as  shown  in 
Fig.  16.  The  transducers  were  mounted  on  a  custom  fabricated  white  ABS  spherical  interface 
with  a  distance  of  71.6mm  from  the  center  of  the  array,  and  formed  into  seven  layers.  Each 
transducer  had  a  central  frequency  of  5MHz  and  a  reception  bandwidth  of  greater  than  80%, 
resulting  in  an  almost  isotropic  spatial  resolution  about  0.2mm.  The  active  area  of  the 
transducer  was  3mm  in  diameter,  and  the  angular  acceptance  was  about  15  degree.  The  signal 
from  the  192  transducers  was  amplified  by  16  homemade  preamplifier  boards  and  coupled  into 
a  64  channel  parallel  data  acquisition  system  with  3:1  multiplexing,  as  indicated  with  three 
different  colors. 

Rats  were  elevated  to  the  center  of  the  spherical  interface  through  a  chamber  fixed  at  the  tank 
bottom,  whose  top  was  about  15mm  beneath  the  interface  center,  and  a  transparent  plastic 
wrap  was  used  to  cover  the  chamber  top.  Ultrasound  gel  was  used  to  couple  the  sound  from  the 
animal  to  the  plastic  wrap.  The  ultrasonic  transducer  array  with  the  spherical  interface,  placed  in 
the  water  bath,  was  therefore  acoustically  coupled  to  the  mouse  tissues.  The  laser  light  was 
delivered  through  a  concave  lens  giving  a  homogenous  illumination  on  the  head  of  the  rat.  In 
this  study,  two  wavelengths  were  used,  1064nm  and  710nm.  The  1064nm  light  was  from  a 
pulsed  Nd:YAG  (5~20ns  pulse  width),  and  the  710nm  light  was  from  a  Q-switched  Nd:YAG 
pumped  Ti:sapphire  laser.  The  final  laser  intensity  on  the  rat  head  was  about  50mJ/cm^  for 
1064nm,  and  4mJ/cm^  for  710nm,  all  below  the  ANSI  exposure  limits.  Both  the  two  lasers 
operate  at  10Hz,  so  that  the  system  could  record  one  complete  set  of  3D  PAT  data  in  0.33s. 
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Two  injection  needles,  inserted  into  the  front  cortex  of  rat  brain  about  3mm  below  the  scalp, 
were  served  as  EEG  electrodes  as  shown  in  Fig.  16(b).  One  electrode  was  picked  for  reference 
and  ground  signal,  while  the  other  electrode  was  to  record  EEG  signal.  The  EEG  signal  was 
amplified  (RA16PA,  TDtucker-Davis  Tech.)  and  recorded  (RZ5  Bioamp  Processor,  TDtucker- 
Davis  Tech.)  with  a  sampling  rate  about  50kHz.  Both  the  two  electrodes  were  firmly  glued  on 
the  rat  head  and  the  wires  were  tightly  attached  to  the  animal  holder  using  glue  tapes,  so  as  to 
avoid  noise  induced  by  any  possible  unnecessary  movement  of  the  rat  or  experiment  operators. 
The  PAT  system  gave  an  output  synchronized  signal  to  the  EEG  system,  so  that  the  two 
systems  were  well  synchronized. 


Fig.  16  Experimental  setup  of  the  combined  system,  (a)  Schematic  of  the  photoacoustic  system;  (b)  Rat  with  two 
electrodes  in  the  brain;  and  the  iaser  iilumination  is  aiso  indicated 

6.  Animal  preparation 

Two  groups  of  rats  (n=16  rats,  each  about  35~45g,  with  8  rats  for  1064nm  and  the  other  8  rats 
for  710nm)  were  used.  The  rats  were  anaesthetized  with  about  0.25ml  Urethane  (0.25g/mL)  and 
got  hair  on  the  head  shaven.  Before  mounting  the  rats  on  the  homemade  plastic  holder,  a 
syringe  containing  0.04g/ml  PTZ  with  a  scalp  needle  was  inserted  in  the  abdomen  and  fixed 
with  glue.  PTZ  was  used  to  induce  generalized  seizure  in  rat  brain,  which  has  already  been 
established.  Then  the  two  electrodes  were  inserted,  and  the  rats  were  elevated  into  the 
transducer  array  center  of  the  photoacoustic  system  for  scanning.  Each  experiment  lasted 
25min,  and  PTZ  was  administrated  with  a  dose  of  0.5mL  at  5min.  All  rats  were  kept  alive 
through  the  experiment.  All  the  animal  procedures  performed  were  in  accordance  with  the 
approved  University  of  Florida  lACUC  protocols. 

C.  Data  processing 

With  a  frame  rate  of  3.3f/s,  5000  complete  sets  of  PAT  data  were  recorded.  Before  doing  3D 
reconstruction  and  processing,  2D  PAT  images  were  reconstructed  with  data  collected  by  the 
64  transducers  in  the  lowest  layer  of  the  array,  which  were  indicated  with  blue  in  Fig.  16.  The  2D 
images  were  carefully  checked  to  avoid  the  possible  motion  artifact  induced  by  the  movement  of 
the  rat  head  during  the  experiment.  3D  images  were  calculated  using  delay&sum  algorithm, 
which  was  accelerated  by  GPU  parallel  technique.  Then  the  negative  pixel  values  in  were  set  to 
zero.  Each  set  of  the  reconstructed  3D  images  had  a  voxel  number  of  201x201x101 
representing  a  10x10x5mm^  volume,  which  could  be  calculated  in  0.97s.  Amira  (from  Visage 
Imaging,  Inc.)  was  used  to  render  the  3D  images. 

There  were  two  sources  of  noise  in  the  EEG  signal,  a  strong  sharp  10Hz  noise  from  the  Q- 
switch  of  the  laser,  and  the  60Hz  alternating  current  noise.  For  de-noising,  the  deteriorated  data 
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by  the  Q-switch  was  replaced  with  adjacent  data,  and  the  60Hz  alternating  current  noise  was 
removed  by  filtering  out  58-62Hz  in  frequency  domain  (Fig.  17).  Then  the  EEG  data  was  applied 
with  a  low  pass  filter  of  100Hz,  and  transformed  into  a  300000  point  signal  for  the  25min  EEG 
recording,  corresponding  to  sampling  rate  of  200Hz. 

For  neurovascular  coupling  study,  we  first  extracted  the  averaged  PAT  signals  of  the  superior 
sagittal  sinus  (SSS)  of  each  rat  from  the  reconstructed  3D  images  and  compared  it  with  the 
corresponding  EEG  signal.  We  choose  the  PAT  signals  of  the  SSS  because  it  has  the  highest 
SNR  and  is  highly  distinguishable  in  the  reconstructed  3D  image.  Because  1064nm  and  710nm 
are  mainly  absorbed  by  oxy-hemoglobin  (Hb02)  and  deoxy-hemoglobin  (HbR),  the  acquired 
results  for  these  two  wavelengths  can  be  taken  as  the  changes  of  Hb02  and  HbR  respectively. 
All  the  acquired  results  were  examined  by  experts  in  the  neurology. 


A.  Data  preprocessing 


Representative  time  frames  of  2D  images  through  the  experiment  are  displayed  in  Fig.  17(a) 
with  a  time  interval  of  Imin,  compared  with  a  rat  brain  photo  taken  after  the  experiment  (Fig. 
17(b)).  Although  the  image  quality  is  not  excellent  due  to  the  limited  transducer  positions  (only 
64  transducers  are  averaged  once),  the  SSS  and  two  transverse  sinuses  (TS)  are  clearly 
reconstructed,  and  some  other  small  blood  vessels  are  also  revealed.  In  Fig.  17(a)  it  is 
observed  that  there  is  no  palpable  shift  of  the  rat  brain  which  is  a  prerequisite  for  the  following 
data  processing.  In  our  experiments,  3  out  of  8  rats  with  1064nm  moved,  but  all  the  rats  in  the 
second  group  remained  still.  We  also  notice  that  it’s  hard  to  tell  when  the  seizure  begins  just  by 
eyes  from  these  2D  images,  so  quantitative  analyze  is  needed. 

EEG  signals  were  recorded  to  confirm  the  onset,  duration,  and  morphology  of  the  seizure,  which 
is  shown  in  Fig.  18(a).  The  EEG  signal  is  stable  before  the  drug  injection  at  5min,  when  isolated 
high  amplitude  spikes  can  be  easily  seen  in  the  EEG  signal.  After  this  spontaneous  fast 
discharges  begin  to  show  and  increase  in  amplitude  in  the  following  2  to  3  minutes.  The  seizure 
onset  comes  about  142  ±  45s  after  the  drug  injection,  which  is  defined  to  be  the  first  signal  that 
exceeds  6  times  the  standard  deviation  (SD)  of  the  baseline  (the  average  signal  in  the  first 
4min).  As  the  seizure  evolves,  the  firing  frequency  generally  gets  lower,  but  the  EEG  signal 
begins  to  decrease  in  and  keeps  firing  through  the  experiment.  Fig.  18(b)  shows  the  EEG  signal 
in  a  short  time  window  (blue),  compared  with  the  original  data  before  preprocessing  (red).  We 
can  see  that  both  the  noises  from  the  Q-switch  and  the  alternating  power  supply  are  removed. 


Fig.  17  Representative  time  frames  of  2D  PAT  images  of  the  rat  brain,  compared  with  an  anatomicai  image,  (a) 
Representative  time  frames  of  2D  PAT  images  of  the  rat  brain  with  1064nm.  The  time  intervai  is  1min,  and  each 
image  covers  lOxIOmm^  with  201x201  pixels,  (b)  Photo  of  rat  brain  with  hair  and  scalp  being  removed.  The  SSS,  TS, 
and  the  reconstruction  area  is  indicated.  Scale  bar  is  5mm  through  the  figure. 
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Fig.  18  EEG  signal  for  PTZ-induced  seizure,  (a)  Typical  EEG  signal  after  processing,  (b)  The  original  (blue)  and 
processed  (red)  EEG  signal  in  a  short  time  wiondow.. 

S.  3D  image  reconstruction 

Fig.  19  shows  the  reconstructed  3D  images  from  different  views.  The  reconstructed  volume  is 
10  X  10  X  5mm^  with  a  0.05mm  voxel  size.  Similar  to  the  2D  images,  the  SSS  and  TS  can  be 
clearly  revealed  both  in  the  3D  images  of  Hb02  and  HbR,  although  the  SNR  of  HbR  is  not  as 
good  as  that  of  Hb02  due  to  relatively  low  laser  intensity  of  710nm,  as  well  as  the  low 
concentration  of  HbR  in  the  blood  vessels  compared  to  that  of  Hb02.  It  is  also  noticed  that  in  the 
3D  images  of  Hb02,  we  can  see  the  Ostia  end  of  the  straight  sinus  (SS)  below  the  SSS  which 
can  not  be  seen  in  the  2D  images.  All  these  demonstrate  the  excellent  positioning  ability  of  our 
3D  PAT  system. 


Fig.  19  Reconstructed  3D  PAT  images  of  the  rat  brain  with  1064nm  from  different  views.  SSS,  TS,  and  SS  can  be 
seen  in  the  image.  The  image  domain  is  lOxlQxSmm^  with  201x201x101  pixels.  Colarbar  is  indicated. 
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C.  Neurovascular  coupling  study 

For  each  set  of  PAT  data,  the  average  signal  in  a  ROI  that  could  enclose  the  SSS  was  extracted 
and  taken  as  the  change  of  the  SSS  (Fig.  20(a)  and  Fig.  20(b)),  which  was  further  converted 
into  percentage  changes  to  the  baseline.  For  all  the  rat  datas  analyzed  (n=14),  the  EEG  signals 
looked  the  same  in  shape  and  timing  (Fig.  20(c)  and  Fig.  20(d)).  This  means  our  PTZ  induced 
generalized  seizure  modal  is  quite  stable  and  the  drug  administration  is  well  controlled,  so  that 
all  the  rat  results  are  comparable.  The  two  PAT  curves  in  Fig.  20(c)  and  Fig.  20(d)  are  also 
representatives.  The  PAT  signals  for  HbR  is  smoothed  using  a  moving  window  of  3  seconds  to 
improve  the  SNR.  The  neurovascular  coupling  of  our  rat  model  can  be  obtained  by  comparing 
these  two  groups  of  PAT  curves  and  the  corresponding  EEG  signals. 


Fig.  20  Typical  results  for  Hb02  (a  and  c)  and  HbR  (b  and  d)  respectively,  with  the  ROI  for  the  SSS  indicated  in  the  3D  PAT  images. 
The  reconstruction  domain  for  the  3D  images  (a  and  c)  are  10x10x5mm^  with  201x201x101  pixels,  c  and  d  are  the  EEG  and  PAT 
signals  of  SSS  for  a  and  b  respectively.  The  EEG  signals  (blue  line)  are  in  arbitrary  units,  and  the  PAT  curves  (red  line)  are 
represented  as  the  change  with  respect  to  the  baseline.  The  whole  seizure  process  after  the  injection  of  PTZ  can  be  divided  into 
three  stages.  The  PAT  signal  for  HbR  has  been  smoothed  using  a  moving  window  of  3  seconds.  The  time  units  in  c  and  d  is  minutes. 

Both  the  EEG  and  PAT  signals  in  these  two  datas  remain  quite  stable  before  the  drug  injection. 
After  that  the  time  course  of  the  experiment  can  be  generally  divided  into  three  stages,  as 
indicated.  In  the  first  stage,  the  EEG  signal  increases  in  amplitude  as  seizure  initiates.  At  the 
same  time,  there’s  a  prominent  sudden  decrease  in  the  Hb02  signal  while  the  HbR  signal  just 
slightly  increases.  This  indicates  that  for  this  period  of  time,  the  increase  of  blood  supply  was 
inadequate  for  the  metabolism  of  neurons.  In  the  second  step,  we  can  see  that  the  EEG  signals 
continue  to  increase  as  the  seizure  evolves.  However,  the  Hb02  signal  begins  to  increase  and 
the  HbR  begins  to  decrease  at  the  same  time.  Thus  we  conclude  that  the  blood  supply  is  greatly 
increased  due  to  dilation  of  arteries  and  increasing  of  blood  flow.  In  the  third  step,  although  the 
EEG  signal  gradually  decreases  in  firing  frequency,  the  Hb02  signal  is  still  decreased,  while  the 
HbR  signal  increases.  This  means  the  blood  flow  is  decreased. 

Through  the  whole  seizure  process,  we  see  the  coincidence  that  the  two  PAT  signals  are  well 
negatively  correlated.  When  the  Hb02  reaches  its  peak  the  HbR  just  drops  to  its  minimum,  or 
otherwise.  To  make  it  more  convincing,  we  calculate  the  average  changes  of  Hb02  and  HbR 
across  all  the  rat  experiments,  as  shown  in  Fig.  21.  The  two  curves  there  are  achieved  by  first 
normalizing  the  data  in  a  single  experiment  to  the  baseline,  and  then  averaging  the  data  of  all 
rats  in  the  same  group.  The  final  data  are  further  averaged  every  SOsecs,  with  the 
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corresponding  SD  indicated.  The  correlation  coefficient  of  these  two  signals  is  calculated  to  be  - 
0.705. 


Fig.  21  Full  timecourses  of  mean  percent  changes  of  the  Hb02  (a)  and  HbR  (b)  signals,  which  are  negatively 
correlated  with  each  other.  The  SD  for  each  signal  is  indicated  every  30s.  Vertical  bar  represents  the  amplitude  of  the 
PAT  signal  changes  to  the  baseline. 

Neurovascular  coupling  study  of  epilepsy  requires  simultaneous  recording  of  neuronal  activity 
and  macroscopic  hemodynamic  signals  in  the  brain,  bringing  significant  challenges  in  terms  of 
experimental  paradigm  and  instrument  performance.  We  successfully  integrated  the  EEG 
system  into  our  real-time  3D  PAT  system  that  is  based  on  a  64-channel  parallel  data  acquisition 
system  and  a  spherical  transducer  array,  and  realized  the  simultaneous  recording  of  neuron 
activity  and  hemodynamic  signal  in  the  whole  rat  brain  over  the  entire  epileptic  seizure  process. 
The  neurovascular  coupling  of  the  intraperitoneal  PTZ  rat  model  here  is  studied  with  the  Hb02 
and  HbR  signals  of  the  SSS.  PTZ  has  been  used  experimentally  to  induce  generalized  seizure, 
where  hemodynamic  changes  can  be  found  all  over  the  brain.  As  one  of  the  main  blood  vessels 
in  the  brain,  SSS  is  highly  involved  in  the  PTZ-induced  seizures.  When  the  seizure  begins,  we 
observed  a  sudden  decrease  in  HbOa  for  1  to  2  minutes,  followed  by  a  dramatic  increase.  The 
‘dip’  observed  here  is  much  similar  to  the  reported  ‘epileptic  dip’  with  ORIS,  but  found  on  the 
cerebral  cortex.  It’s  different  from  the  ‘initial  dip’  which  only  lasts  several  seconds.  In  the  second 
stage  of  the  seizure,  the  Hb02  signal  keeps  on  increasing  and  finally  succeeds  the  baseline, 
which  supports  the  idea  that  the  increase  of  cerebral  blood  flow  oversupplies  the  increased 
metabolic  demands  of  the  neurons  in  epileptic  seizures.  And  the  negative  correlation  we 
observed  between  the  Hb02  and  HbR  signals  in  SSS  here  can  be  explained  by  the  fact  that  the 
veins  don’t  usually  dilate  as  much  as  the  arteries. 

There  are  limitations  in  our  work.  We  only  studied  the  PAT  signals  of  SSS,  but  in-depth 
information  of  the  vascular  coupling  is  more  embedded  in  the  local  hemodynamic  changes. 
Thus  the  signals  of  different  regions  in  the  brain  tissue  need  to  be  extracted,  based  on  which, 
the  neural  network  analysis  can  be  carried  out. 

8.  Miniature  wearabie  PAT  array  system 

We  developed  a  miniature  wearable  PAT  system  that  employs  a  homemade  64-element 
transducer  array. 
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Fig.  22  Miniature  wearable  PAT  imaging  system  and  validation.  (A)  Experimental  photograph  of 
PAT  imaging  system  mounted  on  a  rat  head.  (B)  Corresponding  diagram  showing  the  device 
components.  (C)  Schematic  diagram  of  the  fabrication  of  PVDF  transducer  array.  (D)  Noninvasive 
PAT  image  of  a  freely  moving  rat  brain  with  skin  and  skull  intact.  (D)  Open-skull  photograph  of  the 
rat  cortex  surface  after  the  PAT  experiments 


The  PAT  detecting  module,  shown  in  Fig.22A,  is  composed  with  a  cap-like  housing,  64- 
element  transducer  array,  64  signal  wires  and  4  fiber  bundles.  Fig.22B  shows  the  schematic 
diagram  of  the  cap-like  housing  which  is  designed  and  fabricated  with  3D  printing  technique. 
There  are  four  parts  in  the  housing  structure,  component  1  serves  as  the  bottom  holder  used  to 
support  the  whole  detecting  module  and  also  as  mounting  part  with  the  help  of  two  screws 
placed  oppositely  on  both  side;  component  2  is  used  to  hold  the  transducer  array;  component  3 
is  the  cover  to  protect  the  PAT  imaging  space;  component  4  is  used  to  connect  the  fiber  bundles 
to  the  imaging  area.  As  shown  in  Fig.  22C,  we  utilize  PVDF  film  to  fabricate  the  transducer  array 
using  its  unique  feature  of  flexibility.  In  order  to  make  a  64-element  array  based  on  a 
94.2mmx5mm  size  PVDF  film  (1.47mm  width  and  5mm  height  for  each  element),  we  leave  one 
face  of  the  PVDF  as  a  sharing  ground  and  divided  the  other  face  into  64  rectangles  and  each 
one  is  working  as  an  independent  acoustic  detector.  The  benefit  of  this  design  is  that  the  array 
is  more  electronic  stable  and  it’s  easy  to  be  assembled  into  the  transducer  holder  in  the  housing. 
Fig.22D  shows  a  typical  non-invasive  PAT  image  obtained  with  the  miniature  PAT  imaging 
system  of  a  freely  moving  rat  brain  with  skin  and  skull  intact.  Fig.  22E  is  the  photograph  of  the 
rat  cortex  surface  with  skin  and  skull  removed  after  PAT  experiment.  We  also  built  a  64-channel 
amplify  and  data  acquisition  system  which  can  acquire  64-channel  signals  simultaneously  at  a 
sampling  rate  of  50MHz.  The  imaging  rate  of  the  PAT  system  can  achieve  10  frames/s  without 
averaging. 

8.1  Experiment  and  animal  preparation 

The  bottom  holder  (component  1 )  was  mounted  on  the  rat  head  for  one  hour  before  performing 
experiment  in  order  to  let  the  rat  used  to  the  device.  The  detecting  space  was  filled  with 
phantom  for  acoustic  signal  transports  to  the  transducer  array.  Anti-load  methods  were 
employed  to  reduce  the  weight  on  the  rat  and  it  ensures  the  least  restriction  on  experimental 
animal  goal  of  this  work.  Besides,  EEG  recordings  were  utilized  in  the  experiment  as  a  standard 
reference  to  the  PAT  results.  We  conduct  experiments  on  awake,  freely  moving  rats  (80-1  OOg) 
in  both  resting  and  seizure  (drug  induced  generalized  model.  Pentylenetetrazol,  PTZ,  O.lg/Kg) 
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state  using  1064  nm  wavelength  pulsed  laser  light.  All  procedures  were  approved  by  the 
University  of  Florida  Animal  Care  and  Use  Committee  and  conducted  in  accordance  with  the 
National  Institutes  of  Health  Guide  for  the  Care  and  Use  of  Experimental  Animals. 

8.2  Results 


Fig.  23  shows  the  PAT  and  EEG  results  of  a  behaving  rat  in  resting  state.  The  imaging  time  is  4 
min  and  the  imaging  rate  is  3.33  frames/s  (3  times  averaging).  Fig.  23A  is  the  recorded  EEG 
results  and  Fig.  23B  shows  reconstructed  PAT  signal  change  of  MCV  (middle  cerebral  vessel) 
marked  with  white  arrows  in  Fig.  23C.  Bottom  row  of  Fig.  23C  is  the  close  up  view  of  the  top  row 
PAT  images  at  different  time.  From  Fig.  23A  and  23B  we  can  see  that  there  is  a  corresponding 
significant  change  at  the  time  around  50s  in  both  the  EEG  and  PAT  results.  It  indicates  an 
intense  neural  activity  which  has  been  verified  by  a  moving  of  the  rat  according  to  the  recorded 
video. 


Fig.  23  Noninvasive  PAT  imaging  of  a  freely  moving  rat  in  resting  state.  (A).  EEG  recording;  (B). 
Reconstructed  PAT  signal  changes  of  MCV  marked  with  white  arrows  in  (C);  (C).  PAT  images  of 
the  rat  brain  at  different  time. 


Fig.  24  shows  the  PAT  and  EEG  results  of  an  awake,  freely  moving  rat  in  seizure  state  with 
PTZ  drug  injected  at  the  time  of  90  s.  The  imaging  time  is  20  min  and  the  imaging  rate  is  3.33 
frames/s  (3  times  averaging).  Fig.  24A  is  the  recorded  EEG  results  and  Fig.  24B  shows 
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reconstructed  PAT  signal  changes  of  three  locations  marked  with  white  arrows  in  Fig.  24C.  Fig. 
24C  shows  PAT  images  at  different  time. 
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Fig.  24  Noninvasive  PAT  imaging  of  a  freely  moving  rat  with  drug  induced  seizure.  (A).  EEG  recording;  (B). 
Reconstructed  PAT  signal  changes  MCV  marked  with  white  arrows  in  (C);  (C).  PAT  images  of  the  rat  brain 
at  different  time. 


Fig.  25  shows  the  PAT  and  EEG  results  of  an  anesthetized  rat  in  seizure  state  with  PTZ 
drug  injected  at  the  time  of  180  s.  The  imaging  time  is  20  min  and  the  imaging  rate  is  3.33 
frames/s  (3  times  averaging).  Fig.  25A  is  the  recorded  EEG  results  and  Fig.  25B  shows 
reconstructed  PAT  signal  change  of  MCV  marked  with  white  arrow  in  Fig.  24C.  Fig.  24C 
shows  PAT  images  at  different  time. 
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Fig.  25  Noninvasive  PAT  imaging  of  a  freely  moving  rat  with  drug  induced  seizure.  (A).  EEG  recording;  (B). 
Reconstructed  PAT  signal  changes  of  MCV  marked  with  white  arrows  in  (C);  (C).  PAT  images  of  the  rat 
brain  at  different  time. 


8.3  Discussion 

For  the  first  time,  we  have  presented  a  miniature  wearable  photoacoustic  tomography  system 
for  non-invasive  real-time  imaging  in  freely  moving  rats.  The  system  is  based  on  a  ring-shaped 
PVDF  transducer  array  containing  64  elements.  With  the  64-channel  amplify  and  data 
acquisition  system,  it  can  achieve  imaging  rate  of  10  f/s  which  is  the  frequency  of  the  pulsed 
laser  we  utilized  in  the  experiments.  Both  the  rest  and  drug-induced  seizure  results  demonstrate 
the  imaging  ability  of  our  miniature  PAT  system  on  freely-moving  rats,  as  well  as  the 
neurovascular  effect. 

However,  there  are  still  some  limitations  of  our  current  miniature  PAT  system.  The  first  is  its 
relatively  low  spatial  resolution  due  to  the  limited  sensitivity  of  PVDF  and  the  limit  number  of 
transducer  element.  The  signals  we  can  acquire  are  mainly  from  the  middle  cerebral  blood 
vessel.  The  second  is  the  light  illumination.  We  currently  use  four  fiber  bundles  placed 
according  to  the  shape  of  the  main  blood  vessels  in  the  brain  which  results  in  non-uniform  light 
distribution  on  the  detecting  area.  The  third  but  not  the  last  one,  happened  during  experiments, 
is  the  fixation  of  the  whole  probe.  We  applied  two  screws  held  against  the  skull  of  the  rat  to 
mount  the  probe  on  the  head,  but  it’s  always  not  enough  to  bear  the  rat’s  intensive  motion, 
especially  at  the  seizure  onset  state. 
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Key  Research  Accomplishments 

1.  We  have  completed  the  design  for  an  array  based  PAT  system,  which  is  the  most  important 
part  of  Task  1. 

2.  We  have  developed  three  novel  schemes  that  can  enhance  our  current  reconstruction 
software. 

3.  We  have  completed  the  construction  of  the  proposed  an  array  based  real-time  PAT  system, 
and  calibrated  and  tested  the  system  using  extensive  phantom  experiments. 

4.  We  have  built  the  animal  interface  and  successfully  tested  it  for  in  vivo  imaging  of  rat  brain. 

5.  We  have  developed  three  novel  schemes  that  can  enhance  our  current 
reconstruction  software. 

6.  We  have  conducted  phantom  experiments  that  confirmed  our  software  enhancement. 

7.  We  have  conducted  phantom  experiments  that  confirmed  our  software  enhancement. 

8.  We  have  developed  the  advanced  boundary  conditions  scheme,  and  implemented  the  3D 
reconstruction  codes  and  their  parallelization  which  allowed  efficient  high  performance  PAT 
image  reconstruction. 

9.  We  have  conducted  simulation  and  phantom  experiments  that  confirmed  our  software 
enhancement. 

10.  We  have  performed  extensive  in  vivo  experiments  using  the  rat  epilepsy  model  to  evaluate 
the  PAT  system  we  developed  in  this  project.  The  results  show  that  we  were  able  to  image 
the  functional  anatomy  of  epileptic  foci  noninvasively  and  in  real  time  at  a  spatial 
resolution  of  ~1 50pm. 

11.  We  have  developed  and  validated  a  rat  model  of  temporal  lobe 
epilepsy. 

12.  We  have  performed  extensive  in  vivo  experiments  to  achieve  concurrent  EEG  and  PAT 
recordings.  Interesting  results  are  obtained  that  may  be  significant  for  seizure  prediction  and 
for  neurovascular  coupling  studies. 

13.  We  have  designed  and  constructed  a  light-weighted  animal/transducer  interface  that  can  be 
used  for  chronic  monitoring  of  freely  moving  animals. 

14.  We  analyzed  the  images  from  the  in  vivo  experiments. 

15.  We  built  a  simultaneous  real-time  3D  photoacoustic  tomography  and  EEG  system  for 
neurovascular  coupling  study  in  an  animal  model  of  epilepsy 

16.  We  developed  a  miniature  wearable  PAT  system  that  employs  a  homemade  64-element 
transducer  array 

Reportable  outcomes  (see  the  Appendix  to  this  Summary  Report) 


Conclusions 

We  have  fulfilled  the  statement  of  work  proposed.  The  real-time  PAT  system  completed  in 
Years  1  and  2  provides  us  a  platform  for  performing  the  extensive  phantom  and  animal 
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experiments  proposed  for  the  coming  years.  In  year  3,  we  implemented  reconstruction  codes 
and  we  performed  simulation  and  phantom  experiments  to  evaluate  the  software  developments. 
In  year  4,  the  real-time  PAT  system  constructed  in  this  project  coupled  with  the  light-weighted 
cap  interface  has  provided  us  a  platform  for  performing  extensive  animal  experiments  in  freely 
moving  rats.  In  Year  5  we  completed  the  proposed  Tasks  6  and  7  for  continuous  recording  of 
EEG  and  PAT  over  a  long  period  of  time  to  monitor  seizure  activity  continuously  to  characterize 
changes  during  interictal,  ictal  and  post-ictal  periods  when  a  spontaneous  seizure  occurs.  We 
analyzed  the  images  from  the  in  vivo  experiments  and  write  final  reports.  We  completed  any 
unfinished  experimentation  during  these  months  such  that  we  have  a  sufficient  number  or  rats 
monitored  with  EEG  and  PAT  to  reach  statistical  significance. 
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A  total  variation  minimization  (TVM)-based  finite  element  reconstruction  algorithm  for  photoacoustic 
(PA)  tomography  is  described  in  this  paper.  This  algorithm  is  used  to  enhance  the  quality  of  reconstructed 
PA  images  with  time-domain  data.  Simulations  are  conducted  where  different  contrast  levels  between 
the  target  and  the  background,  different  noise  levels,  and  different  sizes  and  shapes  of  the  target  are 
studied  in  a  30  mm  diameter  circular  heterogeneous  background.  These  simulated  results  show  that 
the  quality  of  the  reconstructed  images  can  be  improved  significantly  due  to  the  decreased  sensitivity 
to  noise  effect  when  the  TVM  is  included  in  the  reconstruction  algorithm.  The  enhancement  is  further 
confirmed  using  experimental  data  obtained  from  several  phantom  experiments  and  an  in  vivo  animal 
experiment.  ©  2011  Optical  Society  of  America 
OCIS  codes:  100.2980,  170.3010,  170.5120,  170.6960. 


1.  Introduction 

Photoacoustic  tomography  (PAT)  is  an  emerging 
noninvasive  imaging  technique  that  combines  the 
merits  of  high  optical  contrast  and  high  ultrasound 
resolution  in  a  single  modality  [1-3].  In  PAT,  the 
Helmholtz-like  photoacoustic  (PA)  wave  equation 
has  been  commonly  used  as  an  accurate  model  for  de¬ 
scribing  laser-induced  acoustic  wave  propagation  in 
tissue.  While  analytical  reconstruction  methods  have 
been  used  for  PA  image  reconstructions,  the  finite 
element  method  (FEM)-based  approach  appears  to 
be  particularly  powerful  in  this  regard  [4,5] .  The  ad¬ 
vantages  of  the  FEM-based  PAT  method  include 
(1)  quantitative  imaging  capability  by  recovering  op¬ 
tical  absorption  coefficient,  (2)  elimination  of  the  as¬ 
sumption  of  homogeneous  acoustic  medium  needed 
in  anal3d:ical  methods,  (3)  accommodation  of  object 
boundary  irregularity,  and  (4)  appropriate  boundary 
conditions  implementations. 

Our  finite  element  PAT  approaches  are  based  on  a 
regularized  iterative  Newton  method,  which  have 
been  tested  and  evaluated  using  extensive  simula- 
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tions  and  phantom  experiments  in  both  the  fre¬ 
quency  and  time  domain  [4-7].  Compared  with  the 
frequency  domain  method,  the  time  domain  method 
can  reconstruct  images  with  less  background  arti¬ 
facts,  especially  when  the  target  size  is  very  small, 
and  can  provide  a  more  accurately  recovered  target 
size  [7].  While  these  results  are  encouraging  and  pro¬ 
mising,  we  realize  that  measurement  noises  are  still 
the  major  factor  affecting  the  quantitative  accuracy 
of  the  reconstructed  images.  For  example,  errors  for 
quantitatively  recovering  optical  absorption  coeffi¬ 
cient  can  be  as  large  as  10%  for  simulated  data  with 
5%  noise  added  and  20%  for  experimental  data  [6,7]. 

In  an  effort  to  reduce  the  noise  effect  and  further 
enhance  the  quantitative  accuracy  of  PA  image 
reconstruction,  in  this  paper,  we  consider  a  total 
variation  minimization  (TVM)  scheme  within  our 
FEM-based  reconstruction  framework.  Our  existing 
reconstruction  algorithms  are  based  on  the  least- 
squares  criteria  (i.e.,  the  regularized  Newton  meth¬ 
od)  [8,9]  that  stand  on  the  statistical  argument  that 
the  least-squares  estimation  is  the  best  estimator 
over  an  entire  ensemble  of  all  possible  pictures.  Total 
variation,  on  the  other  hand,  measures  the  oscilla¬ 
tions  of  a  given  function  and  does  not  unduly  punish 
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discontinuities  Hence,  one  can  hypothesize 

that  a  hybrid  of  these  two  minimization  schemes 
should  be  able  to  provide  higher-quality  image  re¬ 
constructions.  In  fact,  the  strategy  of  finding  minimal 
total  variation  has  proved  to  be  successful  in  applica¬ 
tions  including  electrical-impedance  tomography 
[11],  microwave  imaging  [12],  image  processing 
[W,l^,  W],  optimal  design  [15],  and  diffuse  optical  to¬ 
mography  [1^].  The  TVM  has  also  been  implemented 
and  tested  in  a  non-FEM-based  PA  reconstruction 
framework  [17-19]. 

In  Section  2,  we  describe  in  detail  the  implementa¬ 
tion  of  TVM  within  our  existing  PAT  reconstruction 
framework  in  the  time  domain.  In  Section  3,  we  de¬ 
monstrate  the  enhanced  reconstruction  algorithm 


using  simulated,  phantom,  and  in  vivo  data.  The 
conclusions  are  made  in  Section  4 

2.  TVM  Scheme 

To  describe  our  TVM  method,  we  first  briefiy  intro¬ 
duce  the  FEM-based  regularized  Newton  method  in 
the  time  domain.  The  time  domain  PA  wave  equation 
in  tissue  can  be  described  as  follows  [20]: 

\d‘^p{r,t)  ^{r)pdJ{t) 

C,  3(  ’ 

where  p  is  the  pressure  wave,  Vq  is  the  speed  of  acous¬ 
tic  wave  in  the  medium,  /?  is  the  thermal  expansion 
coefficient,  Cp  is  the  specific  heat,  O  is  the  absorbed 


Fig.  1.  (Color  online)  Reconstructed  absorbed  energy  density  images  from  simulated  data  with  and  without  the  TVM  enhancement  under 
different  noise  levels  (case  1).  (a)  Without  the  TVM,  0%  noise,  (b)  with  the  TVM,  0%  noise,  (c)  without  the  TVM,  10%  noise,  (d)  with  the 
TVM,  10%  noise,  (e)  without  the  TVM,  25%  noise,  (f)  with  the  TVM,  25%  noise.  The  axes  (left  and  bottom)  illustrate  the  spatial  scale  in 
millimeters,  whereas  the  color  scale  (right)  records  the  absorbed  energy  density  in  millijoules  per  cubed  millimeter. 
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energy  density  and  J{t)  =  8{t-tQ)  is  assumed  in 
our  study. 

Expanding  p  as  the  sum  of  coefficients  multiplied 
by  a  set  of  basis  function  y/j,  p  =  Y1  PjWj^  fhe  finite 


Fig.  2.  (Color  online)  Comparison  of  the  exact  and  reconstructed 
absorbed  energy  density  profiles  along  transect  y  =  0  mm  for  the 
images  appearing  in  Fig.  1.  (a)  0%  noise,  (b)  10%  noise,  (c)  25% 
noise. 


element  discretization  of  Eq.  (p  can  then  be  written 
as  [7] 


i/AjVp  •  hAl  = 


i 


dJ 


y/idS. 


(2) 


The  first-order  absorbing  boundary  conditions 
used  are  [21] 


Vp  •  h  =  - 


1  dp 

Vq  dt 


2r’ 


(3) 


where  h  is  the  unit  normal  vector. 

In  both  the  forward  and  inverse  calculations,  the 
unknown  coefficient  O  needs  to  be  expanded  in  a 
similar  fashion  to  p  as  a  sum  of  unknown  parameters 
multiplied  by  a  set  of  locally  spatially  varying 
Lagrange  polynomial  basis  functions.  Thus,  the  ma¬ 
trix  form  of  Eq.  (2)  becomes 


[K]{p}  +  [C]{p}  +  [M]{p}  =  {B},  (4) 


where  the  elements  of  matrix  [K],  [C],  and  [M]  are 


and  the  column  vectors  {p},  {p},  {p},  and  {5}  are 


{p}  =  {pi,P2,---pnV; 
{p}  =  {pi,P2,"-pjv}^; 
{p}  =  {Pi,P2,---PnV- 


Here,  the  Newmark’s  time-stepping  scheme  has  been 
used  for  the  discretization  of  time  dimension  [22,23], 
which  is  a  commonly  used  implicit  method  for  the 
second-order  propagation  equations  such  as  Eq.  (4). 

To  form  an  image  from  a  presumably  uniform  in¬ 
itial  guess  of  the  absorbed  energy  density  distribu¬ 
tion,  we  need  a  method  of  updating  O  from  its 
starting  value.  This  update  is  accomplished  through 
the  least-squares  minimization  of  the  following 
functional: 


M 

P{p,O)  =  ^{p0-pj)2,  (5) 

7=1 

where  p^  and  pj  are  observed  and  computed  acoustic 
field  data  for  j  =  1, 2,  •  •  •  M  boundary  locations.  Using 
the  regularized  Newton  method,  we  obtained  the 
following  equation  for  updating  O: 
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Fig.  3.  (Color  online)  Reconstructed  absorbed  energy  density  images  from  simulated  data  with  and  without  the  TVM  enhancement  with 
different  target  sizes  (case  2).  (a)  Without  the  TVM,  2  mm  diameter  target;  (b)  with  the  TVM,  2  mm  diameter  target;  (c)  absorbed  energy 
density  profiles  along  the  transect  y  =  0  mm  for  the  images  appearing  in  Figs.  1(a)  and  1(b)  (4  mm  diameter  target);  (d)  absorbed  energy 
density  property  profiles  along  the  transect  y  =  0  mm  for  the  images  appearing  in  Figs.  3(a)  and  3(b)  (2  mm  diameter  target).  In  (a)  and  (b), 
the  axes  (left  and  bottom)  illustrate  the  spatial  scale  in  millimeters,  whereas  the  color  scale  (right)  records  the  absorbed  energy  density  in 
millijoules  per  cubed  millimeter. 


(3'^3'  +  /lI)Aj  =  -p'"),  (6) 

where  =  {p\,p%  ■  ■  -phY,  p'^  =  ' '  'PmY’ 

is  the  update  vector  for  the  absorbed  optical  energy 
density  is  the  Jacobian  matrix  formed  by  3p/30  at 
the  boundary  measurement  sites,  X  is  the  regulariza¬ 
tion  parameter  determined  by  combined  Marquardt 
and  Tikhonov  regularization  schemes,  and  I  is  the 
identity  matrix. 

Two  typical  approaches  exist  for  minimizing  total 
variation:  a  constrained  minimization  through  the 
solution  of  the  nonlinear  PA  equation  [8,^]  and  an 
unconstrained  minimization  by  addition  of  the  total 
variation  as  a  penalty  term  to  the  least-squares  func¬ 
tional  [W,  1^,14,^].  From  a  computational  stand¬ 
point,  unconstrained  minimizations  are  much 
easier  to  implement  and  require  fewer  modifications 
to  the  existing  algorithm  [13].  Thus,  in  this  study,  the 
unconstrained  TVM  is  used. 

We  now  incorporate  the  total  variation  of  O  as  a 
penalty  term  by  defining  a  new  functional  [9,  W,^] : 

F(p,0)=F(p,0)+L(0).  (7) 


Here, 


L(0)  =  J  +  (5^drdy 


(8) 


is  the  penalty  term,  and  and  S  are  typically 
positive  parameters  that  need  to  be  determined  nu¬ 
merically.  The  minimization  of  Eq.  (7)  proceeds  in 
standard  fashion  by  the  differentiation  of  F  with  re¬ 
spect  to  each  nodal  parameter  that  constitutes  the  O 
distribution;  simultaneously,  all  these  relations  are 
set  to  zero.  These  steps  lead  to  the  following  non¬ 
linear  system  of  equations 


dF  ^ 

-pY  ^ 


where 


(9) 


X  dxdy. 
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Fig.  4.  (Color  online)  Reconstructed  absorbed  energy  density  images  from  simulated  data  with  and  without  the  TVM  enhancement  with 
different  contrast  levels  between  the  target  and  the  background  (case  3).  (a)  Without  the  TVM,  1.5 : 1  contrast;  (b)  with  the  TVM,  1.5 : 1 
contrast;  (c)  absorbed  energy  density  profiles  along  the  transect  y  =  0  mm  for  the  images  appearing  in  Figs.  1(a)  and  1(b)  (2 : 1  contrast); 
(d)  absorbed  energy  density  profiles  along  the  transect  y  =  0  mm  for  the  images  appearing  in  Figs.  4(a)  and  4(b)  (1.5 : 1  contrast).  In  Figs.  4 
(a)  and  4(b),  the  axes  (left  and  bottom)  illustrate  the  spatial  scale  in  millimeters,  whereas  the  color  scale  (right)  records  the  absorbed  energy 
density  in  millijoules  per  cubed  millimeter. 


Similar  to  Eq.  (6),  the  following  matrix  equation  for 
TVM  constrained  inversion  can  be  obtained: 

+  i?  +  i/)  A/  =  - pc)  _  (10) 


where 


30]^  3^2  30jY 
.  ^^2 
30]^  3O2  30jY 


30i  3O2  '  '  ’  30^ 


and  V  =  {Vi,V2,---V;vr. 


3.  Results  and  Discussion 

In  this  section,  the  TVM-enhanced  reconstruction 
algorithm  is  tested  and  evaluated  using  both 
simulated  and  experimental  data.  For  comparative 
purposes,  reconstructions  without  the  TVM  enhance¬ 
ment  are  also  presented. 


A.  Simulations 

In  these  simulations,  a  dual-meshing  method,  as  de¬ 
tailed  elsewhere  [6] ,  is  used  for  fast  yet  accurate  in¬ 
verse  computation.  The  fine  mesh  used  for  the 
forward  calculation  consisted  of  3627  nodes  and 
7072  elements,  while  the  coarse  mesh  used  for  the 
inverse  calculation  had  930  nodes  and  1768  ele¬ 
ments.  All  the  images  obtained  from  the  method 
without  the  TVM  are  the  results  of  three  iterations, 
while  those  obtained  from  the  TVM-enhanced  meth¬ 
od  are  the  results  of  15  or  more  iterations.  Parallel 
code  was  used  to  perform  these  calculations  on 
Beowulf  clusters  with  8  central  processing  units 
(CPUs).  The  computational  time  can  be  further  re¬ 
duced  if  clusters  with  more  than  8  CPUs  are  used.  As 
mentioned  earlier,  the  parameters  and  5  were 
determined  through  numerical  experimentation.  A 
constant  value  of  ^  =  0.001  was  sufficient  for  the  cur¬ 
rent  simulation  and  experimental  studies,  while  the 
value  of  0)^  is  related  to  the  signal-to-noise  ratio  of 
the  measurements  [W].  For  the  simulations  pre¬ 
sented,  co^  =  0.5  for  cases  1,  2,  and  3,  and  = 
1.0  for  case  4  was  used. 

For  the  first  simulation,  the  test  geometry  is  a 
30  mm  diameter  circular  background  containing  an 
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Fig.  5.  (Color  online)  Reconstructed  absorbed  energy  density 
images  from  simulated  data  with  and  without  the  TVM  enhance¬ 
ment  for  three  targets  having  different  shapes  (case  4).  (a)  Exact 
image,  (b)  without  the  TVM,  (c)  with  the  TVM,  (d)  absorbed  energy 
density  profiles  along  the  transect  y  =  0  mm.  In  (a)-(c),  the  axes 
(left  and  bottom)  illustrate  the  spatial  scale  in  millimeters, 
whereas  the  color  scale  (right)  records  the  absorbed  energy  density 
in  millijoules  per  cubed  millimeter. 


off-center  4  mm  diameter  target  region.  The  target 
had  O  =  2.0mJ/mm^,  while  the  background  had 
O  =  l.OmJ/mm^.  In  this  case,  the  image  reconstruc¬ 
tion  was  performed  with  0%,  10%,  and  25%  noise 
added  to  the  ''measured”  data.  Figure  1  gives  two  sets 
of  absorbed  energy  density  images  recovered  using 
the  reconstruction  method  without  [Figs.  1(a),  1(c), 
and  1(e)]  and  with  the  TVM  enhancement  [Figs.  1(b), 
1(d),  and  1(f)]  under  the  conditions  of  different  noise  le¬ 
vels.  As  can  be  seen,  enhancement  of  the  reconstruc¬ 
tion  by  the  incorporation  of  the  TVM  is  obvious  over 
the  method  without  the  TVM.  To  provide  a  more  quan¬ 
titative  assessment  of  these  images.  Fig.  2  is  included, 
in  which  the  reconstructed  absorbed  energy  density 
profiles  are  displayed  along  transects  through  the 
target  for  the  images  shown  in  Fig.  1.  We  find  that 
the  TVM-enhanced  method  not  only  can  improve 
the  quality  of  the  recovered  images,  but  also  can  de¬ 
crease  the  sensitivity  of  the  method  to  noise  effect. 

The  second  simulation  is  intended  to  investigate 
how  the  target  size  affects  the  TVM  enhancement. 
In  this  case,  no  noise  was  added  to  the  "measured” 
data,  and  the  diameter  of  the  off-center  target  was 
2  mm.  Figures  3(a)  and  3(b)  give  two  sets  of  absorbed 
energy  density  images  recovered  using  the  method 
without  the  TVM  [Fig.  3(a)]  and  with  the  enhance¬ 
ment  [Fig.  3(b)],  while  Figs.  3(c)  and  3(d)  present 
the  quantitative  profiles  of  the  absorbed  energy  den¬ 
sity  along  transects  through  the  target  for  the  images 
shown  in  Figs.  1(a),  1(b),  3(a),  and  3(b).  Again,  con¬ 
siderable  improvement  can  be  observed  from  the 
reconstructed  results  when  the  TVM  is  invoked  com¬ 
pared  with  the  method  without  the  TVM.  It  is  also 
interesting  to  note  that  the  enhancement  for  this 
case  is  more  striking  than  that  for  the  first  one  where 
the  background  contained  a  larger  target. 

The  third  simulation  aims  to  see  how  the  contrast 
level  of  the  absorbed  energy  density  between  the  target 
and  the  background  impacts  the  TVM  enhancement. 
In  this  case,  no  noise  was  added  to  the  "measured” 
data,  and  the  off-center  target  had  a  diameter  of 
4  mm  and  0  =  1.5  m J / mm^ .  Figures  4(a)  and  4(b)  pre¬ 
sent  two  sets  of  absorbed  energy  density  images  recov¬ 
ered  using  the  method  without  the  TVM  [Fig.  4(a)]  and 
with  the  TVM  enhancement  [Fig.  4(b)] ,  while  Figs.  4(c) 
and  4(d)  show  the  quantitative  profiles  of  the  absorbed 
energy  density  along  transects  through  the  target  for 
the  images  shown  in  Figs.  1(a),  1(b),  4(a),  and  4(b).  We 
can  see  that  the  images  formed  by  incorporation  of  the 
TVM  are  clearly  enhanced  qualitatively  in  visual  con¬ 
tent  relative  to  that  obtained  using  the  method  with¬ 
out  the  TVM.  We  also  note  that  lower  the  contrast  level 
is,  the  larger  the  enhancement. 

In  the  fourth  simulation,  three  targets  having  dif¬ 
ferent  shapes  (circular,  elliptical,  and  rectangular) 
were  embedded  in  the  background,  and  the  absorbed 
energy  density  of  these  targets  was  2.5,  1.5,  and 
2.0mJ/mm^,  respectively.  A  noise  level  of  25%  was 
added  to  the  "measured”  data  in  this  case.  Figure  5 
shows  the  exact  and  the  recovered  absorbed  energy 
density  images  as  well  as  the  quantitative  absorbed 
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Fig.  6.  UQI  calculated  from  the  recovered  images  with  and  without  TVM  enhancement  from  simulated  data,  (a)  Case  1  when  different 
noise  levels  are  considered,  (b)  cases  2-4. 


Fig.  7.  (Color  online)  Reconstructed  absorbed  energy  density  images  from  the  three  phantom  experiments,  (a)  Case  1  without  the 
TVM,  (b)  case  1  with  the  TVM,  (c)  case  2  without  the  TVM,  (d)  case  2  with  the  TVM,  (e)  case  3  without  the  TVM,  (f)  case  3  with  the 
TVM.  The  axes  (left  and  bottom)  illustrate  the  spatial  scale  in  millimeters,  whereas  the  color  scale  (right)  records  the  absorbed  energy 
density  in  millijoules  per  cubed  millimeter. 
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energy  density  property  profiles  along  the  transect 
that  across  these  targets.  Again,  the  improvement 
in  image  quality  is  apparent. 

In  order  to  quantitatively  evaluate  the  reconstruc¬ 
tion  quality  using  the  reconstruction  method  with 
and  without  TVM  enhancement,  we  used  the  univer¬ 
sal  quality  index  (UQI)  to  measure  the  degree  of  si¬ 
milarity  between  the  reconstructed  and  exact  images 
[26].  We  first  interpreted  an  image  f  as  vectors  of  size 
N  :  f  =  (/i,  Ar  •  where  superscript  T  denotes 

the  matrix  transpose,  and  N  denotes  the  number 
of  image  data  acquired  from  the  FEM-based  algo¬ 
rithm.  We  then  defined  image  means,  variances  and 
covariances  over  the  whole  image  domain  as 

(11) 

k=l 


= 


N 


k=i 


(12) 


where  j  =  0  and  1,  and 


Cov{fi,fO}  =  -p).  (13) 

^  k=l 

Hence,  the  UQI  is  defined  as 

„  2Cov(fM<')  2/‘y 

UQI  measures  the  image  similarity  between  the  re¬ 
constructed  (f^)  and  reference/exact  (f^)  images,  and 
its  value  ranges  between  0  and  1.  The  value  of  the 
UQI  is  closer  to  1  when  the  reconstructed  image  is 
more  similar  to  the  exact  image.  We  calculated  the 
UQIs  for  the  simulation  cases  presented  earlier, 
and  the  results  are  given  in  Fig.  6,  where  Figs.  6(a) 
and  6(b)  show  the  UQIs  for  case  1  when  different  noise 
level  is  considered  and  for  all  the  four  cases  without 
added  noise.  In  Fig.  6,  the  horizon  axis  indicates  the 
case  number  (from  cases  1  to  4,  presented  earlier),  and 


X  position  (mm) 


Fig.  8.  (Color  online)  Recovered  absorbed  energy  density  profiles  along  (a)  y  =  -7.0nini  crossing  the  3  mm  diameter  target  for 
experimental  case  1,  (b)  y  =  8.0  mm  crossing  the  2  mm  diameter  target  for  experimental  case  1,  (c)  y  =  1.0  mm  for  experimental  case 
2,  and  (d)y  =  6.5  mm  for  experimental  case  3. 
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Fig.  9.  (Color  online)  CNR  calculated  for  the  recovered  images 
using  the  method  with  and  without  TVM.  (a)-(c)  Images  shown 
in  Figs.  7(b),  7(d),  and  7(f)  with  the  ROIs  marked  (four  pairs  of 
ROIs:  solid  line  circle,  ^-ROI;  dashed  line  circle,  6-ROI).  (d)  CNR 
computed  for  the  four  pairs  of  ROIs  shown  in  (a)-(c).  In  Figs.  9 
(a)-9(c),  the  axes  (left  and  bottom)  illustrate  the  spatial  scale  in 
millimeters,  whereas  the  color  scale  (right)  records  the  absorbed 
energy  density  in  millijoules  per  cubed  millimeter. 


the  vertical  axis  shows  the  value  of  UQI.  We  immedi¬ 
ately  note  from  Fig.  6  that  the  TVM-based  method 
provides  significantly  better  UQIs  than  that  without 
TVM. 

B.  Phantom  and  In  Vivo  Experiments 

In  this  section,  both  phantom  and  in  vivo  experimen¬ 
tal  data  were  used  to  confirm  our  findings  from  the 
simulations.  The  experimental  setup  used  for  collect¬ 
ing  the  phantom  data  was  a  pulsed  NdiYAG  laser- 
based  single  transducer  (IMHz)  scanning  system, 
which  was  described  in  detail  elsewhere  [5].  Three 
phantom  experiments  were  conducted.  In  the  first 
two  experiments,  we  embedded  one  or  two  objects 
with  a  size  ranging  from  3  to  0.5  mm  in  a  50  mm  dia¬ 
meter  solid  cylindrical  phantom.  The  phantom  mate¬ 
rials  used  consisted  of  Intralipid  as  a  scatterer  and 
India  ink  as  an  absorber  with  Agar  powder  (1%- 
2%)  for  solidifying  the  Intralipid  and  India  ink  solu¬ 
tion.  The  absorption  coefficient  of  the  background 
phantom  was  0.01  mm“^,  while  the  absorption 
coefficient  of  the  target(s)  was  0.03  mm“^.  In  the  last 
experiment,  we  used  a  single  target-containing 
phantom,  aiming  to  test  the  capability  of  detecting 
targets  having  low  optical  contrasts  relative  to  the 
background  phantom.  In  this  case,  the  target  had 
an  absorption  coefficient  of  0.015  mm“^.  The  reduced 
scattering  coefficients  of  the  background  phantom 
and  targets  were  1.0  and  3.0  mm“^  for  the  first  two 
experiments,  and  1.0  and  2.0  mm“^  for  the  last 
experiment. 

A  total  of  120  receivers  were  equally  distributed 
along  the  surface  of  the  circular  background  region. 
In  the  reconstructions,  the  fine  mesh  used  for  the  for¬ 
ward  calculation  consisted  of  5977  nodes  and  11,712 
elements,  while  the  coarse  mesh  used  for  the  inverse 
calculation  had  1525  nodes  and  2928  elements.  The 
reconstructed  images  were  the  results  of  three  and 
15  iterations  for  the  method  without  and  with  the 
TVM,  respectively.  For  these  experimental  cases, 
(D^  =  1.0  and  5  =  0.001  for  the  cases  1  and  3  and 
0)^  =  2.0  and  5  =  0.001  for  case  2  were  used. 

Here  we  note  that,  in  the  experimental  situation, 
the  effect  of  ultrasonic  transducer  response  should  be 
considered  in  the  PAT  image  reconstruction  because 
the  transducer  mechanical-electrical  impulse  re¬ 
sponse  as  well  as  the  transducer  aperture  effect 
may  introduce  significant  errors  in  the  forward  mod¬ 
el.  To  reduce  these  effects  in  our  calculations,  we  ap¬ 
plied  the  normalized  acoustic  pressure  distribution 
in  the  reconstruction  and  used  the  appropriate 
parameters  (initial  values,  etc.)  obtained  from  an 
optimization/searching  method  [27]. 

Figure  7  shows  the  reconstructed  absorbed  energy 
density  images  for  all  the  three  experimental  cases, 
while  Fig.  8  presents  quantitative  absorption  coeffi¬ 
cient  profiles  along  transects  through  one  target  for 
the  images  shown  in  Fig.  7.  We  see  that  considerably 
enhanced  images  are  achieved  with  the  TVM,  espe¬ 
cially  when  the  target  is  small  (case  2),  or  when  the 
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(a)  (b) 

Fig.  10.  (Color  online)  Recovered  absorbed  energy  density  images  from  the  rat  brain  (a)  without  the  TVM,  (b)  with  the  TVM.  The  axes  (left 
and  bottom)  illustrate  the  spatial  scale  in  millimeters,  whereas  the  color  scale  (right)  records  the  absorbed  energy  density  in  millijoules  per 
cubed  millimeter. 


contrast  level  between  the  target  and  the  back¬ 
ground  is  low  (case  3). 

Because  there  is  no  true  or  exact  image  available 
for  the  experimental  cases,  we  used  the  contrast-to- 
noise  ratio  (CNR)  to  quantitatively  evaluate  the  ex¬ 
perimental  results.  To  compute  CNR,  we  selected 
two  regions  of  interest  (ROIs),  each  having  the  same 
size.  One  region  is  in  the  target  (referred  as  ^-ROI) 
and  the  other  is  in  the  background  (denoted  as  b- 
ROI).  The  CNR  is  defined  as 


CNR  = 


l^(^)   f{b)  I 


(15) 


where  and/^^^  respectively,  denote  the  mean  over 
the  ^-ROI  and  6-ROI,  and  denotes  the  variances 
over  the  6-ROI.  We  used  solid  and  dashed  circles  to 
express  the  ^-ROI  and  6-ROI  for  the  three  experimen¬ 
tal  cases.  We  calculated  the  CNR  for  each  pair  of  the 
t-  and  6-ROIs  using  Eq.  (15)  and  present  the  results 
in  Fig.  9(d).  We  see  that  the  CNR  obtained  using  the 
method  with  TVM  is  clearly  larger  than  that  using 
the  method  without  TVM,  indicating  that  the  TVM- 
based  method  is  less  sensitive  to  the  noise  effects. 

As  a  final  example,  our  methods  were  tested  using 
in  vivo  data  collected  previously  from  an  animal  (rat) 
model  of  epilepsy  [28].  Focal  seizures  were  induced 
by  microinjection  of  bicuculline  methiodide  (BMI) 
into  the  parietal  neocortex.  The  in  vivo  data  were  col¬ 
lected  using  the  same  pulsed  NdiYAG  laser-based 
scanning  PAT  system.  In  the  reconstruction,  the  fine 
mesh  used  for  the  forward  calculation  consisted  of 
17,713  nodes  and  34,944  elements,  while  the  coarse 
mesh  used  for  the  inverse  calculation  had  4489  nodes 
and  8736  elements.  The  in  vivo  images  obtained  were 
the  results  of  two  and  five  iterations  for  the  methods 
without  and  with  the  TVM,  respectively.  We  found 
that  =  0.1  and  ^  =  0.001  appeared  to  provide 
optimal  results. 

Figure  W  gives  the  recovered  images  for  the  rat 
brain  scanned  10  min  after  the  injection  of  BMI  with¬ 
out  and  with  the  use  of  the  TVM.  The  in  vivo  results 
shown  here  confirm  that  both  of  the  reconstruction 


methods  can  provide  quality  images  (the  arrow  indi¬ 
cates  the  detected  seizure  focus),  while  notable  en¬ 
hancement  in  image  quality  can  be  observed  with 
the  TVM-based  method.  For  example,  the  actual  con¬ 
tinuous  middle  cerebral  artery  (indicated  by  the 
dashed  circle)  is  shown  discontinued  in  Fig.  10(a) 
(without  the  TVM),  while  this  feature  is  clearly 
depicted  in  Fig.  10(b)  (with  the  TVM). 

4.  Conclusions 

We  have  presented  a  time-domain  FEM-based  PA 
image  reconstruction  method  that  incorporates  the 
TVM.  The  results  shown  in  this  work  indicate  that 
the  TVM-based  method  is  able  to  offer  clear  enhance¬ 
ment  in  image  reconstruction  over  the  method  with¬ 
out  the  TVM,  not  only  in  terms  of  the  location,  size, 
and  shape  of  the  target,  but  also  in  terms  of  the  ab¬ 
sorbed  energy  density  property/optical  absorption 
coefficient  values  themselves.  In  addition,  the  results 
have  shown  that  the  inclusion  of  the  TVM  in  our  re¬ 
construction  algorithm  is  highly  effective  in  the  pre¬ 
sence  of  noisy  data.  Finally,  it  is  important  to  note 
that  the  TVM-based  method  may  be  ideal  for  image 
reconstruction  involving  a  low  contrast  level  between 
the  target  and  the  background  or  small  targets  em¬ 
bedded  in  the  background. 

This  research  was  supported  in  part  by  the  Depart¬ 
ment  of  Defense  Congressionally  Directed  Medical 
Program. 
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Abstract:  Photoacoustic  tomography  (PAT)  is  an  emerging  non-invasive 
imaging  technique  with  great  potential  for  a  wide  range  of  biomedical 
imaging  applications.  However,  the  conventional  PAT  reconstruction 
algorithms  often  provide  distorted  images  with  strong  artifacts  in  cases 
when  the  signals  are  collected  from  few  measurements  or  over  an  aperture 
that  does  not  enclose  the  object.  In  this  work,  we  present  a  total-variation- 
minimization  (TVM)  enhanced  iterative  reconstruction  algorithm  that  can 
provide  excellent  photoacoustic  image  reconstruction  from  few-detector  and 
limited-angle  data.  The  enhancement  is  confirmed  and  evaluated  using 
several  phantom  experiments. 
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OCXS  codes:  (100.2980)  Image  enhancement;  (170.3010)  Image  reconstruction  techniques; 
(170.5120)  Photoacoustic  imaging;  (170.6960)  Tomography 
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1.  Introduction 

Photoaeoustie  tomography  (PAT)  is  an  emerging  non-invasive  imaging  teehnique  that 
eombines  the  merits  of  high  optieal  eontrast  and  high  ultrasound  resolution  in  a  single 
modality  [1-3].  In  PAT,  the  Helmholtz-like  photoaeoustie  wave  equation  has  been  eommonly 
used  as  an  aeeurate  model  for  deseribing  laser-indueed  aeoustie  wave  propagation  in  tissue. 
While  a  variety  of  analytie  PAT  reeonstruetion  algorithms  have  been  developed  [4-6],  the 
finite  element  method  (FEM)  based  approaeh  appears  to  be  partieularly  powerful  in  this 
regard,  beeause  it  ean  provide  quantitative  imaging  eapability  by  reeovering  optieal  absorption 
eoeffieient  [7,8],  eliminate  the  assumption  of  homogeneous  aeoustie  medium  needed  in 
analytieal  methods,  aeeommodate  objeet  boundary  irregularity  and  allow  appropriate 
boundary  eonditions  implementations. 

A  praetieal  need  exists  for  reeonstruetion  of  photoaeoustie  images  from  few 
measurements,  as  this  ean  greatly  reduee  the  required  seanning  time  and  the  number  of 
ultrasound  sensors  plaeed  near  or  on  the  boundary  of  an  objeet  to  reeeive  the  laser-indueed 
aeoustie  signals.  In  addition,  in  many  praetieal  implementations  of  PAT  the  photoaeoustie 
signals  are  reeorded  over  an  aperture  that  does  not  enelose  the  objeet,  whieh  results  in  a 
limited-angle  tomographie  reeonstruetion  problem.  In  sueh  eases,  the  existing  reeonstruetion 
algorithms,  whieh  are  based  on  the  least-squares  eriteria  (i.e.,  the  regularized  Newton  method) 
[9,10],  often  generate  distorted  images  with  severe  artifaets.  Total-variation-minimization 
(TVM),  on  the  other  hand,  has  proved  to  be  a  powerful  tool  for  limited-data  image 
reconstruction  in  diffraction  tomography  and  computed  tomography  (CT)  [11,12].  TVM  based 
PAT  algorithms  have  also  been  implemented  and  tested  using  numerically  simulated  limited- 
view  data  [13-16].  Here  we  describe  an  unconstrained  TVM  FEM-based  iterative 
reconstruction  algorithm,  and  for  the  first  time  present  experimental  evidence  that  the  TVM- 
based  algorithm  offers  excellent  photoaeoustie  image  reconstruction  from  few-detector  and 
limited-angle  data. 

This  paper  is  organized  as  follows.  In  Section  2,  the  FEM  based  PAT  reconstruction 
algorithm  and  the  unconstrained  TVM  scheme  are  reviewed  briefly.  The  experimental 
validation  of  our  TVM  FEM-based  algorithm  is  presented  in  Section  3.  The  conclusions  are 
made  in  Section  4. 

2.  Method 

We  first  briefly  describe  the  existing  FEM-based  photoaeoustie  reconstruction  algorithm 
detailed  elsewhere  [17].  The  time-domain  photoaeoustie  wave  equation  in  tissue  can  be 
described  as  follows: 


Be  Bt  ’ 

where  p  is  the  pressure  wave;  Vq  is  the  speed  of  acoustic  wave  in  the  medium;  P  is  the  thermal 
expansion  coefficient;  is  the  specific  heat;  O  is  the  absorbed  energy  density; 
j(t)  =  5 is  assumed  in  our  study. 
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To  form  an  image  from  a  presumably  uniform  initial  guess  of  the  absorbed  energy  density 
distribution,  we  need  a  method  of  updating  O,  the  desired  quantity  from  its  starting  value.  This 
update  is  accomplished  through  the  least-squares  minimization  of  the  following  functional: 

M  . 

(2) 

y=i 

where  and  are  the  measured  and  computed  acoustic  field  data  for  /  =  1,2,---M 

boundary  locations.  Using  the  regularized  Newton  method,  we  obtained  the  following  matrix 
equation  for  updating  O: 

(3"3  +  Al)A;r  =  3"(p‘’-;>^),  (3) 

where  P°  =iyP° P°m^  and  p‘'  =[P[,P2,'"Pm)  is  the  update  vector  for  the 

absorbed  optical  energy  density;  3  is  the  Jacobian  matrix  formed  by  at  the  boundary 

measurement  sites;  X  is  the  regularization  parameter  determined  by  combined  Marquardt  and 
Tikhonov  regularization  schemes;  and  I  is  the  identity  matrix. 

We  now  incorporate  the  total  variation  of  O  as  a  penalty  term  by  defining  a  new  functional 
[18]: 


F(;?,0)  =  F(;?,0)  +  Z(0).  (4) 

Here  dxdy  is  the  penalty  term,  and  co^  and  S  are  typically  positive 

parameters  that  need  to  be  determined  numerically.  The  minimization  of  Eq.  (4)  can  be 
realized  by  the  differentiation  of  F  with  respect  to  each  nodal  parameter  that  constitutes  the  O 
distribution  and  by  setting  each  of  the  resulting  expression  to  zero,  leading  to  the  following 
system  of  equations: 


y=i 


(5) 


where  V.  -dLjd^b.  .  Again,  through  the  regularized  Newton  method,  the  following  matrix 
equation  for  TVM-based  inversion  is  obtained  [18]: 


(6) 


where  V  is  formed  by  5Z/50  and  R  is  formed  by  5E/50  .  At  this  point,  our  solution 
procedure  for  solving  Eq.  (6)  and  the  regularization  parameter  selection  are  identical  to  those 
described  previously.  Hence,  it  becomes  clear  that  the  only  additions  to  our  new  algorithm 
result  from  the  assembly  of  matrix  R  and  the  construction  of  column  vector  V. 

3.  Results 

In  this  section  our  TVM  enhanced  reconstruction  algorithm  is  tested  and  evaluated  using 
phantom  experimental  data.  For  comparative  purposes,  reconstruction  results  without  the 
TVM  enhancement  are  also  presented. 

The  experimental  setup  used  for  collecting  the  phantom  data  was  a  pulsed  ND:  YAG  laser 
based  single  transducer  (IMHz)  scanning  system,  which  was  described  in  detail  elsewhere  [8]. 
Three  phantom  experiments  were  conducted.  In  the  first  two  experiments,  we  embedded  one 
or  two  objects  with  a  size  ranging  from  3  to  0.5  mm  in  a  50  mm-diameter  solid  cylindrical 
phantom.  The  absorption  coefficient  of  the  background  phantom  was  0.01  mm~\  while  the 
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absorption  coefficient  of  the  target(s)  was  0.03  In  the  last  experiment,  we  used  a  single- 

target-containing  phantom,  aiming  to  test  the  capability  of  detecting  target  having  low  optical 
contrasts  relative  to  the  background  phantom.  In  this  case,  the  target  had  an  absorption 
coefficient  of  0.015  mm~^  The  reduced  scattering  coefficients  of  the  background  phantom  and 
targets  were  1.0  and  3.0  mm“^  for  the  first  two  experiments,  and  1.0  and  2.0  mm”^  for  the  last 
experiment.  In  the  reconstructions,  we  used  a  dual-meshing  scheme  [19]  where  the  fine  mesh 
used  for  the  forward  calculation  consisted  of  5977  nodes  and  11712  elements,  while  the 
coarse  mesh  used  for  the  inverse  calculation  had  1525  nodes  and  2928  elements.  All  the 
images  obtained  from  the  method  without  the  TVM  are  the  results  of  three  iterations,  while 
those  obtained  from  the  TVM-enhanced  method  are  the  results  of  twenty  or  more  iterations. 
Parallel  code  was  used  to  speed  up  these  calculations  on  Beowulf  clusters  with  8  CPUs.  The 
computational  speed  can  be  further  increased  if  clusters  with  more  than  8  CPUs  are  used.  As 
mentioned  above,  the  parameters  and  S  were  determined  through  numerical 

experimentation.  From  our  experience,  a  constant  value  of  (5  =  0.001  was  sufficient  for  the 
current  experimental  studies,  while  the  value  of  is  related  to  the  signal-to-noise  ratio  of  the 
measurements.  For  our  experimental  cases,  =1.0  and  S  =  0.001  for  the  first  case,  and 

=  0.5  and  ^  =  0.001  for  the  second  and  third  cases  were  used. 

The  reconstruction  results  based  on  few-detector  data  from  the  three  experimental  cases 
are  shown  in  Figs.  1,  2  and  3,  respectively.  In  each  figure,  the  images  in  the  top  and  bottom 
rows,  respectively,  present  the  recovered  absorbed  energy  density  images  using  the  algorithm 
without  and  with  the  TVM  where  120,  60,  30,  and  15  detectors  were  equally  distributed  along 
the  surface  of  the  circular  background  region.  As  expected,  the  conventional  algorithm 
without  the  TVM  provides  high  quality  images  only  when  the  number  of  detectors  was 
relatively  sufficient  (Figs.  la,b.  Figs.  2a,b,  and  Figs.  3a,b).  The  images  reconstructed  from 
few-detector  data  by  the  conventional  algorithm  without  the  TVM  contained  severe  artifacts 
and  distortions  (Figs.  lc,d.  Figs.  2c,d,  and  Figs.  3c,d).  Considerably  enhanced  images  are 
achieved  using  the  method  with  the  TVM,  especially  when  the  number  of  the  detectors  is 
reduced  to  30  or  15  (Figs.  lg,h.  Figs.  2g,h,  and  Figs.  3g,h).  We  also  note  that  the 
computational  efficiency  of  our  TVM  based  algorithm  for  recovering  a  small  absorber  (Fig.  2) 
and  a  larger  absorber  (Fig.  3)  is  similar  when  the  number  of  the  detectors  is  insufficient. 

Images  reconstructed  from  the  limited-angle  data  for  the  first  case  using  the  two  methods 
are  displayed  in  Fig.  4  where  the  top  and  bottom  rows,  respectively,  show  the  recovered 
absorbed  energy  density  images  using  the  algorithm  without  and  with  the  TVM  when  120 
detectors  over  360°,  60  detectors  over  180°,  and  30  detectors  over  90°,  were  equally 


Fig.  1.  Reconstructed  photoacoustic  images  based  on  few-detector  data  for  case  1.  (a),  120 
detectors,  without  TVM.  (b),  60  detectors,  without  TVM.  (c),  30  detectors,  without  TVM.  (d), 
15  detectors,  without  TVM.  (e),  120  detectors,  with  TVM.  (f),  60  detectors,  with  TVM.  (g),  30 
detectors,  with  TVM.  (h),  15  detectors,  with  TVM. 
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Fig.  2.  Reconstructed  photoacoustic  images  based  on  few-detector  data  for  case  2.  (a),  120 
detectors,  without  TVM.  (b),  60  detectors,  without  TVM.  (c),  30  detectors,  without  TVM.  (d), 
15  detectors,  without  TVM.  (e),  120  detectors,  with  TVM.  (Q,  60  detectors,  with  TVM.  (g),  30 
detectors,  with  TVM.  (h),  15  detectors,  with  TVM. 


Fig.  3.  Reconstructed  photoacoustic  images  based  on  few-detector  data  for  case  3.  (a),  120 
detectors,  without  TVM.  (b),  60  detectors,  without  TVM.  (c),  30  detectors,  without  TVM.  (d), 
15  detectors,  without  TVM.  (e),  120  detectors,  with  TVM.  (Q,  60  detectors,  with  TVM.  (g),  30 
detectors,  with  TVM.  (h),  15  detectors,  with  TVM. 


Fig.  4.  Reconstructed  photoacoustic  images  based  on  limited-angle  data  for  case  1.  (a),  120 
detectors  over  360°,  without  TVM.  (b),  60  detectors  over  180°,  without  TVM.  (c),  30  detectors 
over  90°,  without  TVM.  (d),  120  detectors  over  360°,  with  TVM.  (e),  60  detectors  over  180°, 
with  TVM.  (f),  30  detectors  over  90°,  with  TVM. 

distributed  along  the  surface  of  the  circular  background  region.  Again,  strong  artifacts  and 
distortions  exist  in  the  images  recovered  with  the  method  without  the  TVM  when  the  data 
collected  is  angle-limited  (Figs.  4b, c),  while  considerably  improved  images  are  reconstructed 
using  the  TVM  enhanced  algorithm  (Figs.  4e,f). 
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In  order  to  evaluate  quantitatively  the  reconstruction  quality  using  the  method  with  and 
without  TVM  enhancement,  we  use  the  following  universal  quality  index  (UQI)  [20]: 


2C<,..{f'.f“}  2/7 


Here  the  image  f  can  be  interpreted  as  vectors  of  size  N:  f  =  (/p/2,---/^)^  ,  where  N  denotes 
the  number  of  image  data  acquired  from  the  FEM  based  algorithm.  is  the  image  means, 
(j^  is  the  variances,  and  C(9v|f\f^|  is  the  covariances  over  the  whole  image  domain,  where  j 

=  0  and  1 .  UQI  measures  the  image  similarity  between  the  reconstructed  (  f  ^ )  and  reference 
(  f  ^ )  images,  and  its  value  ranges  between  0  and  1 .  The  value  of  the  UQI  is  closer  to  1  when 
the  reconstructed  image  is  more  similar  to  the  reference  image. 

We  calculated  the  UQIs  for  the  three  cases  presented  above,  and  the  results  are  given  in 
Fig.  5  where  Figs.  5a  and  5b  show  the  results  from  the  few-detector  data  based  on  the  method 
without  and  with  the  TVM,  respectively,  while  Fig.  5c  presents  the  results  from  the  limited- 
angle  data  for  the  first  case  based  on  the  two  methods.  Here  the  image  data  recovered  by  120 
detectors  was  regarded  as  the  reference  one.  The  computed  UQIs  shown  in  Fig.  5  confirm  the 
observation  that  the  TVM  enhanced  PAT  algorithm  provides  significantly  better  image  quality 
compared  to  the  conventional  method  without  the  TVM. 


(a)  (c) 

Fig.  5.  UQIs  calculated  from  the  reeovered  images  for  the  three  eases  based  on  few-deteetor 
data  without  TVM  (a)  and  with  TVM  (b),  and  for  ease  1  based  on  limited  angle  data  with  and 
without  TVM. 


4.  Conclusions 

In  this  work,  we  have  implemented  and  evaluated  an  unconstrained,  total-variation- 
minimization  method  for  time-domain  FEM  based  PAT.  This  study  has  confirmed  that  in  the 
situation  that  the  data  is  collected  from  few  measurements  or  over  an  aperture  that  does  not 
enclose  the  object,  the  developed  TVM  based  PAT  algorithm  provides  considerably  improved 
image  reconstruction  compared  to  the  conventional  methods.  The  application  of  this  TVM 
enhanced  reconstruction  algorithm  will  significantly  reduce  the  number  of  ultrasound 
transducers  and  scanning  time  needed  for  high  quality  photoacoustic  image  reconstruction. 
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Abstract:  A  real-time  three-dimensional  (3D)  photoacoustic  imaging 
system  was  developed  for  epilepsy  imaging  in  small  animals.  The  system  is 
based  on  a  spherical  array  containing  1 92  transducers  with  a  5  MHz  central 
frequency.  The  signals  from  the  192  transducers  are  amplified  by  16 
homemade  preamplifier  boards  with  26  dB  and  multiplexed  into  a  64 
channel  data  acquisition  system.  It  can  record  a  complete  set  of  3D  data  at  a 
frame  rate  of  3.3  f/s,  and  the  spatial  resolution  is  about  0.2  mm.  Phantom 
experiments  were  conducted  to  demonstrate  the  high  imaging  quality  and 
real  time  imaging  ability  of  the  system.  Finally,  we  tested  the  system  on  an 
acute  epilepsy  rat  model,  and  the  induced  seizure  focus  was  successfully 
detected  using  this  system. 

©2012  Optical  Society  of  America 

OCXS  codes:  (110.5120)  Photoacoustic  imaging;  (170.6920)  Time-resolved  imaging. 
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1.  Introduction 

Photoacoustic  tomography  (PAT)  is  a  hybrid  method  that  is  capable  of  imaging  optical 
absorption  of  tissue  through  the  detection  of  ultrasound  waves  generated  by  a  short  laser  pulse 
due  to  transient  thermoelastic  expansion.  It  has  an  imaging  resolution  that  is  superior  to  pure 
optical  imaging  at  centimeter  scale  depths.  To  date  PAT  has  been  applied  to  the  detection  of 
breast  cancer,  skin  cancer  and  osteoarthritis  in  humans  [1-3],  and  functional  brain  imaging  in 
small  animals  [4,5]. 

For  small  animal  brain  imaging,  most  prior  PAT  studies  were  mostly  based  on  a  single¬ 
transducer  scanning  system  or  a  circular  or  linear  array  of  transducers  for  only  2D  imaging 
[5-7],  and  the  scanning  in  z  direction  is  needed  for  3D  imaging  purpose,  which  leads  to  a  non- 
optimal  elevational  resolution  in  the  z  direction.  2D  planar  array  of  transducers  has  recently 
been  employed  for  3D  PAT  imaging  [8,9].  However,  due  to  the  limited  aperture  of  2D  planar 
array  transducers,  features  with  high  aspect  ratio  or  with  orientations  oblique  to  the  transducer 
surface  suffer  from  distortion,  and  the  azimuthal  resolution  is  reduced.  Thus,  2D  planar  array 
of  transducers  is  not  suitable  for  small  animal  brain  imaging. 

Compared  to  a  2D  planar  array,  a  spherical  array  can  offer  more  complete  angular  views 
of  the  object,  providing  both  high  resolution  and  accurate  feature  definition  regardless  of 
shape  or  location  of  the  object.  The  use  of  sparse  spherical  arrays  for  3D  PAT  imaging  has 
been  recently  reported  [1,10],  but  due  to  the  insufficient  number  of  transducers  used,  these 
arrays  were  not  designed  for  small  animal  brain  imaging.  The  goal  of  this  work  is  to  present  a 
sparse  spherical  array  based  PAT  system  that  is  specifically  designed  for  real  time  3D  imaging 
of  small  animal  brains.  The  work  is  a  natural  extension/improvement  of  our  previous  work 
which  reported  for  the  first  time  2D  PAT  imaging  of  epileptic  focus  in  small  animals  [5].  The 
current  PAT  system  is  based  on  a  2D  spherical  array  of  196  transducers  coupled  with  parallel 
data  acquisition,  offering  a  temporal  resolution  of  0.33  s  for  data  acquisition.  We  demonstrate 
this  system  using  static/dynamic  phantom  and  in  vivo  animal  experiments.  To  the  best  of  our 
knowledge,  this  is  the  first  work  reporting  3D  PAT  imaging  of  epileptic  focus  in  small 
animals. 


2.  System  descriptions 


Figure  1  depicts  the  block  diagram  of  our  real-time  3D  PAT  system.  A  Ti: sapphire  laser 
optically  pumped  with  a  Q-switched  Nd:YAG  laser  sent  8-12  ns  pulses  at  10  Hz  with  a 
wavelength  tunable  from  690  to  1015  nm.  The  beam  was  delivered  with  an  optical  fiber 


Fig.  1 .  Block  diagram  of  our  real-time  3D  PAT  system.  The  inset  is  a  photograph  of  the  elose- 
up  view  of  the  ehamber  holding  the  rat  head. 
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through  an  opening  on  the  top  of  the  transducer  array  and  produced  an  approximately  uniform 
illumination  in  a  2cm-diameter  area  onto  the  sample. 

The  transducer  array  consisted  of  192  transducers  placed  along  a  custom  fabricated  white 
ABS  spherical  interface  containing  610  through  holes  with  counter  bores,  as  shown  in  Fig. 
2(a).  More  holes  were  drilled  so  that  the  selection  of  the  transducer  positions  on  the  ball  can 
be  flexible.  The  interface  has  an  outer  diameter  of  160  mm  and  an  inner  diameter  of  140  mm, 
and  the  diameter  of  the  holes  in  the  ball  is  5.7  mm,  which  fitted  well  with  the  transducer  (5.5 
mm  outer  diameter).  Each  transducer  (Custom  designed  from  Blatek,  Inc.)  has  a  central 
frequency  of  5  MHz  with  a  reception  bandwidth  of  greater  than  80%.  The  active  area  of  the 
transducer  is  3  mm  in  diameter  and  the  angular  acceptance  is  about  15  degree.  The  transducers 
were  glued  onto  the  interface  with  epoxy  which  can  be  removed  to  allow  the  position  change 
of  the  transducers. 

There  were  1 6  preamplifier  boards  separately  sealed  in  4  metal  boxes,  and  each  board  had 
12  input  channels,  4  output  channels,  and  2  digital  signal  control  inputs.  Within  each  board, 
12  dedicated  operational  amplifier  modules  (AD8099)  individually  amplify  the  input  signals 
with  a  fixed  gain  of  26  dB,  and  then  the  amplified  signals  were  multiplexed  into  the  4  output 
channels  by  4  multiplexer  chips  (MAX  4051),  which  were  controlled  by  a  USB  10  digital 
module  (USB-1024LS,  Measurement  Computing)  through  the  2  digital  inputs. 

The  64-channel  parallel  data  acquisition  system  consisted  of  eight  8-channel  PCI  cards 
(PCIAD850,  US  Ultratek)  in  an  industrial  computer.  For  each  channel,  3000  sampling  points 
were  collected  at  50  MHz  sampling  rate  in  10  bits,  and  stored  in  a  32k  on  board  memory 
before  they  were  transferred  to  the  host  machine.  Amplifiers  with  a  programmable  gain  of  0  to 
80  dB  along  with  a  16  MHz  low  band  pass  filter  were  built  into  the  data  acquisition  system.  A 
Labview  program  controlled  the  data  acquisition,  and  the  acquired  data  was  stored  on  hard 
disk  for  further  image  processing.  Images  were  reconstructed  with  a  delay-and-sum  algorithm 
[11].  For  3D  image  display,  the  reconstructed  results  were  normalized  to  0~255  after  setting 
the  negative  values  to  be  zero.  Then  3D  images  were  rendered  with  Amira  (from  Visage 
Imaging,  Inc.)  with  different  thresholds  as  indicated  in  the  colorbar  shown  in  each  image. 


Fig.  2.  The  spherical  transducer  array,  (a)  Photograph  of  the  transducer  array,  (b)  3D  schematic 
of  the  transducer  distribution  on  the  interface. 


The  system  allows  the  selection  of  transducer  positions  on  the  spherical  interface.  The 
total  610  holes  formed  1 1  evenly  spaced  layers  along  the  vertical  direction  of  the  ball,  and  the 
192  transducer  positions  were  indicated  with  three  different  colors  in  Fig.  2b.  Three  different 
colors  were  used  to  indicate  the  3  to  1  multiplexing  from  the  192  transducers  to  the  64- 
channel  data  acquisition.  This  system  can  also  be  used  as  a  real-time  2D  system  operating  at 
10  f/s  using  the  64  transducers  arranged  in  the  vertical  center  layer  (blue).  For  in  vivo 
experiments,  the  rat  head  was  elevated  to  the  center  of  the  spherical  interface  through  a 
chamber  fixed  at  the  tank  bottom,  whose  top  was  about  1 5  mm  beneath  the  interface  center, 
and  a  transparent  plastic  wrap  was  used  to  cover  the  chamber  top.  For  phantom  experiments,  a 
homemade  silicone  holder  was  used  to  hold  the  phantom. 
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3.  Phantom  experiments 

Three  different  types  of  phantoms  were  used:  One  containing  a  point  object  for  system 
calibration,  one  phantom  containing  three  hairs  tilted  along  different  orientations  for  static 
imaging,  and  one  phantom  with  ink  flowing  through  a  thin  tube  embedded  in  a  phantom  for 
real-time  3D  imaging. 

The  point  object  used  for  calibration  was  a  small  spherical  graphite  particle  (0.1  mm  in 
diameter)  located  at  the  center  of  the  spherical  array  and  ensured  an  isotropic  acoustic 
emission  profile  for  all  directions.  We  measured  and  compensated  the  delay  of  time  for  all  the 
192  channels  in  the  radical  direction,  and  reconstructed  this  point  object  after  calibration.  The 
hairs-containing  phantom  was  used  to  demonstrate  the  high  imaging  quality  of  our  system. 
Finally,  we  imaged  an  embedded  tube  filled  with  flowing  ink  to  show  the  real-time  imaging 
ability  of  our  system.  The  tube  had  a  0.3  mm  inner  diameter  and  was  horizontally  placed  in  a 
phantom.  No  averaging  of  signals  was  performed  for  the  phantom  experiments  except  for  the 
point  object  experiment  where  10  times  averaging  was  applied. 

3.1.  System  characterization 

We  calibrated  the  system  by  recording  the  emission  profiles  from  the  spherical  graphite 
particle  for  all  the  192  channels,  and  then  measured  and  compensated  the  time  delay  of  each 
channel.  We  then  evaluated  the  system  resolution  by  reconstructing  the  image  of  the  point 
object. 

200 
100 
0 


Fig.  3.  (a)  and  (c):  x-y  and  z-x  cross  section  images  through  the  eenter  plan  of  the  point  objeet. 
(b)  and  (d):  the  profile  extraeted  in  x  and  z  direetions  from  (a)  and  (e),  respeetively.  Units  are  in 
mm. 


Figures  3a  and  3c  present  the  reconstructed  x-y  and  z-x  cross-section  images  of  the  point 
object  located  at  the  array  center.  The  quality  of  these  images  is  determined  by  both  the 
distribution  and  the  characteristics  of  the  transducers.  The  profiles  of  the  two  reconstructed 
images  were  also  extracted  in  x  and  z  directions,  as  shown  in  Figs.  3b  and  3d,  respectively. 
The  full  width  at  half  maximum  (FWHM)  of  the  profiles  was  measured  to  be  0.19  mm  (x 
direction)  for  Fig.  3b,  and  0.27  mm  (z  direction)  for  Fig.  3d,  compared  to  the  theoretical  value 
of  0.16  mm  for  the  5  MHz  central  frequency  transducer  with  an  estimated  cut  off  frequency  of 
7  MHz.  It  is  noted  that  the  profile  in  Fig.  3d  is  noisier  than  that  in  Fig.  3b.  This  along  with  a 
larger  FWHM  from  Fig.  3d  was  due  to  the  asymmetric  distribution  of  the  transducers.  For 
targets  located  away  from  the  center  of  the  array,  the  radial  resolution  will  stay  nearly  the 
same  as  that  for  a  centrally  located  target,  while  the  lateral  resolution  will  be  linearly  reduced 
with  increased  distance  away  from  the  array  center.  In  our  system,  the  lateral  resolution  will 
be  reduced  by  0.1  mm  when  the  target  is  located  5  mm  off  the  array  center. 

3.2.  Static  phantom  experiments 

Figure  4a  is  the  photograph  of  the  phantom  containing  three  hairs  tilted  along  different 
orientations,  and  Figs.  4b  and  4c  show  the  reconstructed  3D  images  from  two  different  views. 
The  reconstructed  volume  is  10x10x10  mm  with  a  0.1  mm  voxel  size.  The  spatial 
distribution  and  tails  for  all  the  three  hairs  were  clearly  revealed.  This  result  indicates  that  our 
system  is  capable  of  three-dimensionally  imaging  small  objects  of  different  spatial  distribution 
and  orientation  in  high  quality. 
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Fig.  4.  (a):  photograph  of  the  phantom  containing  three  tiled  hairs;  (b)-(c):  reconstructed  3D 
images  of  the  three  hairs  in  two  different  views.  Scale  bar  represents  5  mm. 


3.3.  Dynamic  phantom  experiments 

The  reconstructed  3D  images  of  ink  flowing  through  a  thin  tube  are  shown  in  Fig.  5.  The 
image  domain  is  15x5x10  mm  with  a  0.1  mm  voxel  size.  Figure  5a  is  the  photograph  of  the 
phantom  containing  a  tube  filled  with  ink,  and  Figs.  5b-5j  are  the  reconstructed  3D  images  at 
different  time  points.  The  time  interval  between  two  consecutive  images  was  0.3  s.  These  3D 
images  clearly  tracked  the  flow  through  the  tube  over  the  course  of  2.4  seconds  with  high 
spatial  and  temporal  resolution,  and  the  flowing  speed  of  the  ink  was  measured  to  be  6  ±  0.9 
mm/s  from  the  reconstructed  results  with  a  time  interval  of  0.3  s. 


Fig.  5.  Reconstructed  3D  images  of  ink  flowing  through  a  0.3  mm-tube  embedded  in  a 
background  phantom,  (a):  photograph  of  the  phantom  containing  the  tube,  (b)-(j):  reconstructed 
3D  images  at  different  time  points.  The  time  interval  is  0.3  s.  Scale  bar  represents  5  mm. 


4.  Rat  epilepsy  experiment 

Epilepsy  is  a  serious  brain  disorder  involving  intensive  hemodynamic  changes,  which 
provides  high  endogenous  contrast  for  PAT  imaging  due  to  the  strong  absorption  of  blood  at 
visible  and  NIR  wavelengths.  Compared  with  current  existing  neuron-imaging  methods  (such 
as  MRI,  CT,  PET,  and  SPECT),  PAT  provides  not  only  high  ultrasound  resolution  and  high 
optical  contrast,  but  also  unprecedented  advantage  of  high  temporal  resolution  over  these 
methods,  which  is  critical  for  capturing  seizure  dynamics. 

2D  PAT  of  seizure  focus  on  an  acute  seizure  rat  model  was  demonstrated  for  the  first  time 
by  our  lab  [5],  but  the  observation  of  hemodynamic  changes  during  seizure  onset  was 
hindered  by  the  long  time  scanning  of  a  single  transducer.  Here  we  test  our  real-time  3D 
system  using  the  same  animal  model  to  show  the  hemodynamic  changes  and  reveal  the  3D 
structures  in  the  rat  brain. 

Two  small  rats  (~40g)  were  imaged  with  intact  skull  and  skin  but  hairs  on  the  head  were 
removed.  The  rats  were  anaesthetized  and  mounted  on  the  homemade  plastic  chamber/holder. 
Focal  seizure  was  induced  by  microinjection  of  10  pi  of  1.9  mM  bicuculline  methiodide 
(BMI)  into  the  neocortex  of  one  rat,  while  saline  solution  was  injected  into  the  brain  of 
another  rat  as  control.  In  each  experiment,  the  rat  was  elevated  to  the  transducer  array  center 
and  kept  alive  under  the  water  tank  through  the  whole  experiment.  The  incident  energy  of  the 
730  nm  light  was  maintained  at  8mJ/cm^  below  the  safety  standard.  Seizure  process  was 
recorded  for  50  minutes,  and  the  measurement  from  the  control  rat  was  recorded  for  3  s.  All 
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Fig.  6.  (a)  and  (c):  photograph  of  the  rat  with  BMI  injection  and  the  control  rat  after  scalp 
removed,  (b)  and  (d):  3D  PAT  images  at  6  time  points  for  (a)  and  (c)  respectively.  The  three 
main  blood  vessels  are  indicated  by  the  white  arrows,  and  the  seizure  focus  is  indicated  by  the 
circle  in  (b).  The  time  interval  between  two  successive  images  is  0.3  s.  Scale  bar  represents  10 
mm. 


animal  procedures  were  performed  in  accordance  with  the  approved  University  of  Florida 
lACUC  protocols. 

Figure  6b  presents  the  3D  images  for  the  rat  with  BMI  injection  during  the  seizure  onset  at 
6  time  points,  compared  with  that  for  the  control  rat  in  Fig.  6b.  The  corresponding 
photographs  of  the  two  rats  with  scalp  removed  right  after  the  experiments  are  shown  in  Figs. 
6a  and  6c,  respectively.  The  reconstructed  domain  for  the  images  shown  in  both  Figs.  6b  and 
6d  is  20  X  20  X  4.5  mm,  with  a  0.1  mm  pixel  size.  The  three  main  blood  vessels  on  the  rat 
brain  are  clearly  revealed  for  both  cases,  as  indicated  by  the  white  arrows,  and  for  the  rat  with 
BMI  injection  a  seizure  focus  can  be  clearly  seen  right  at  the  BMI  injection  position  (the  white 
circle  in  Fig.  6b),  where  strong  absorption  is  observed  during  seizure  onset.  The  rapid  changes 
of  the  absorption  both  in  the  main  blood  vessel  and  in  the  seizure  focus  were  observed  from 
Fig.  6b,  while  no  such  changes  were  noted  for  the  control  rat  (Fig.  6c).  The  seizure  focus  had 
a  diameter  of  ~3  mm,  which  is  consistent  with  the  results  published  before  [5].  This 
experiment  demonstrates  that  our  system  can  be  used  to  investigate  the  hemodynamic  changes 
in  small  animal  brain  both  spatially  and  temporally  during  seizure  onset,  although  the  complex 
microvasculature  cannot  be  resolved  due  to  the  limited  number  of  transducers  and  the  simple 
backprojection  reconstruction  method  used  here.  The  micro  vasculature  can  be  revealed  if 
sophisticated  reconstruction  methods  such  as  the  finite  element  based  algorithms  coupled  with 
the  total  variation  minimization  scheme  are  used  [12,13]. 

5  Conclusions 

We  have  presented  a  real-time  photoacoustic  system  for  three-dimensionally  imaging  focal 
cortical  seizures  in  a  rat.  The  system  is  based  on  a  spherical  array  containing  192  discrete 
transducers.  With  the  64-channel  data  acquisition  system  coupled  with  3:1  multiplexing,  it  can 
achieve  a  frame  rate  of  3.3  f/s  with  a  spatial  resolution  of  0.2  mm.  The  3D  imaging 
performance  of  the  system  was  demonstrated  by  both  static  and  dynamic  phantom 
experiments.  We  have  also  tested  our  system  using  an  acute  epilepsy  rat  model  and  obtained 
3D  images  showing  the  hemodynamic  changes  during  seizure  onset. 
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Purpose:  Recently  reported  quantitative  photoacoustic  tomography  (PAT)  has  signihcantly  expanded 
the  utilities  of  PAT  because  it  allows  for  recovery  of  tissue  optical  absorption  coefficient  which  di¬ 
rectly  correlates  with  tissue  physiological  information.  However,  the  recovery  of  optical  absorption 
coefficient  by  the  existing  quantitative  PAT  approaches  strongly  depends  on  the  accuracy  of  absorbed 
energy  density  distribution,  and  on  the  knowledge  of  accurate  strength  and  distribution  of  incident 
light  source.  The  purpose  of  this  study  is  to  develop  a  new  algorithm  for  the  reconstruction  of  optical 
absorption  coefficient  that  does  not  depend  on  these  initial  parameters. 

Methods:  Here  the  authors  propose  a  novel  one-step  reconstruction  approach  that  can  directly  recover 
optical  absorption  coefficient  from  photoacoustic  measurements  along  boundary  domain.  The  authors 
validate  the  method  using  simulation  and  phantom  experiments. 

Results:  The  authors  have  demonstrated  experimental  evidence  that  it  is  possible  to  directly  recover 
optical  absorption  coefficient  maps  using  boundary  photoacoustic  measurements  coupled  with  the 
photon  diffusion  equation  in  just  one  step.  The  authors  found  that  the  method  described  is  able  to 
quantitatively  reconstruct  absorbing  objects  with  different  sizes  and  optical  contrast  levels. 
Conclusions:  Compared  to  the  authors’  previous  two-step  methods,  the  reconstruction  re¬ 
sults  obtained  here  show  that  the  one- step  scheme  can  signihcantly  improve  the  accuracy 
of  absorption  coefficient  recovery.  ©  2012  American  Association  of  Physicists  in  Medicine. 
[http://dx.doi.org/10. 1118/1 .476098 1  ] 
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Biomedical  photoacoustic  tomography  (PAT),  as  a  future 
imaging  modality,  can  visualize  the  internal  structure  and 
function  of  soft  tissues  in  multiscale  (from  millimeters  to  cen¬ 
timeters)  with  high  spatial  resolution  and  excellent  optical 
contrast.  While  conventional  PAT  can  image  tissues  with 
high  spatial  resolution,  it  recovers  only  the  distribution  of  ab¬ 
sorbed  light  energy  density  that  is  the  product  of  tissue  optical 
absorption  coefficient  and  local  optical  huence.  Recently  re¬ 
ported  quantitative  PAT  (qPAT)  (Refs.  4-14)  has  signihcantly 
expanded  the  utilities  of  PAT  as  it  allows  for  recovery  of  tis¬ 
sue  absorption  coefficient  which  directly  correlates  with  tis¬ 
sue  physiological  information.  Ah  these  methods  for  qPAT 
require  a  model-based  reconstruction  algorithm  or  model-free 
method  to  hrst  capture  the  map  of  absorbed  energy  density  or 
initial  pressure.  Then  under  various  assumptions,  the  absorp¬ 
tion  coefficient  is  deduced  from  the  obtained  absorbed  energy 
density/initial  pressure  through  the  photon  diffusion  equation 
or  a  measurement  technique  that  can  provide  the  distribution 
of  exciting  huence. 

However,  there  are  several  fundamental  limitations  associ¬ 
ated  with  the  existing  qPAT  techniques.  First,  one  has  to  know 
the  exact  boundary  rehection  coefficients  as  well  as  the  exact 
strength  and  distribution  of  a  specihed  incident  light  source 
associated  with  the  photon  diffusion/transport  model.  It  re¬ 
quires  careful  experimental  calibration  procedures  in  order  to 


obtain  these  initial  parameters.  Second,  the  recovered  optical 
absorption  coefficient  also  strongly  depends  on  the  accuracy 
of  the  distribution  of  absorbed  energy  density.  Generally,  it  is 
a  challenging  task  to  obtain  an  exact  distribution  of  absorbed 
energy  density  due  to  several  reasons:  the  limited  bandwidth 
of  ultrasound  transducers  for  ultrawideband  photoacoustic 
frequency  response,  dependency  on  the  orientation,  size,  and 
shape  of  the  targets  with  respect  to  the  ultrasound  receiving 
aperture,  and  effect  of  the  optical  scattering.  Finally,  an  op¬ 
tical  measurement  method  for  measuring  the  optical  huence 
will  increase  the  complexity  and  cost  of  hardware  systems. 

To  overcome  the  limitations  mentioned  above,  here  we 
propose  a  novel  one- step  reconstruction  approach  that  can  di¬ 
rectly  recover  optical  absorption  coefficient  from  photoacous¬ 
tic  measurements  along  boundary  domain.  We  validate  the 
method  using  simulations  as  well  as  phantom  experiments. 

The  photoacoustic  wave  equation  in  frequency  domain  for 
an  acoustically  homogeneous  medium  is  written  as  (note  that 
the  method  described  here  can  be  easily  extended  to  time  do¬ 
main  or  acoustically  heterogeneous  cases), 

V  p{r,  0))  -h  k^pir,  co)  =  iko -  or 

Cp 

V'^p{r,oS)  +  klp{r,(o)  =  ik(3^ff-  (1) 

Lp 
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where  p  is  the  pressure  wave;  =  coIcq  is  the  wave  number 
described  by  the  angular  frequency,  co  and  the  speed  of  acous¬ 
tic  wave  in  the  medium,  cq;  is  the  thermal  expansion  coef¬ 
ficient;  Cp  is  the  specific  heat;  E  is  the  absorbed  light  energy 
density,  which  is  the  product  of  absorption  coefficient  /x^  and 
optical  fiuence  O  {E  =  /x^O).  Based  on  the  finite  element  so¬ 
lution  to  Eq.  (1),^  we  use  Marquardt-Tikhonov  regularization 
based  Newton  method  to  directly  recover/update  absorption 
coefficient  from  an  initial  guess  of  /x^^o  by  minimizing  an  ob¬ 
ject  function  composed  of  a  sum  of  the  squared  difference  be¬ 
tween  computed  and  measured  acoustic  pressure  around  the 
boundary  areas  given  by 

L  M 

Min:  F  =  Y,Y.  ^2) 

j  =  l  /  =  1 


And  photon  fiuence  O  in  Eq.  (7)  can  be  computed  by  the  finite 
element  solution  to  the  following  photon  diffusion  equation: 

V  .  Z)(r)VO(r)  -  /x,(r)ch(r)  =  -5(r)  (8) 

in  which  S(r)  is  the  normalized  light  source  term,  D(r)  is  the 
diffusion  coefficient  and  is  considered  a  constant  here.  The 
finite  element  discretized  form  of  Eq.  (8)  is  written  as, 

[K]n.n  .  (9) 

Based  on  Eq.  (9),  the  sensitivity  90/9/Xa  in  Eq.  (7)  is  obtained 
from  the  following  equation  in  the  whole  problem  domain  us¬ 
ing  adjoint  method: 


1— U- 

'  dK' 

1  dp^a  j 

_dlXa. 

As  a  consequence,  for  each  of  the  L  frequencies,  the  complex 
matrix  equation  for  the  forward  problem  in  Eq.  (1)  is  stated 
as 

AnxnP  n  =  (3) 


And  the  following  matrix  equation  capable  of  inverse  so¬ 
lution  of  /Xq,  is  obtained  for  the  inverse  calculation:^ 

(3^3  +  A./)Ax  /tfl  = /tfl.o  +  (4) 

where  p°  =  Y!j=\  PiU)^  •  •  •  ^  PM(j)f  and 

=  T.)=i{p\{j),  P2U),  •  •  • ,  PmU)/  in  which  com- 
plex  p^(j)  and  Pfij)  are  observed  and  computed  normalized 
acoustic  measurements  for  /  =  1 ,  2,  . . . ,  M  boundary  lo¬ 
cations  with  frequency  7,^  and  N  is  the  nodal  number  of 
the  finite  element  mesh  in  the  entire  problem  domain;  A/ 
(dimension:  N)  is  the  update  vector  for  /x^;  ^  (dimension: 
(L  X  2M)  X  AO  is  the  Jacobian  matrix  formed  by  dp! dp. a  at 
the  boundary  measurement  sites;  k  is  a  scalar;  ^  is  calculated 
from  common  backtracking  line  search and  I  is  the  identity 
matrix.  In  consideration  of  £"  =  /x^O,  the  sensitivity  of  pres¬ 
sure  to  absorption  coefficient  (dpfdpia)  equals  the  product 
of  the  sensitivity  of  pressure  to  energy  density  (dpIdE)  in 
full  field  and  the  sensitivity  of  energy  density  to  absorption 
coefficient  (dEldp^a)  in  local  field.  Specifically,  the  elements 
in  Jacobian  matrix  ^  are  determined  by 

dpi 
^  l^aj 


(/  =  1,2,  ...,M;  j,k  =  l,2,  ...,  A) 

(5) 


in  which  the  sensitivities  of  dpIdE  in  each  iteration  can  be 
calculated  from  Eq.  (3)  based  on  the  adjoint  method:^ 


[A] 


\dp] 

~dA~ 

laEj 

laEj 

{P]- 


(6) 


Noted  here  [dA/dE]  =  0  and  the  derivatives  dEldpia  in  Eq.  (5) 
are  further  written  as 

dEk 

^  l^aj 


(^k  +  l^a.kid^k/dl^aj)  = 

\^^a,k(9^k/^^^a.j)  if 


Thus  the  absorption  coefficient  distribution  can  be  recon¬ 
structed  through  a  Newton  iterative  solution  procedure  de¬ 
scribed  by  Eq.  (3)  (forward  solution)  and  Eq.  (4)  (inverse 
solution)  to  minimize  the  objective  function  in  Eq.  (2). 
Equations  (5)-(10)  are  employed  to  calculate  the  Jacobian 
matrix  in  Eq.  (4),  and  Eqs.  (2)-(4)  are  solved  separately  for 
each  iteration.  The  regularization  parameter  is  determined  by 
combined  Marquardt  and  Tikhonov  regularization  schemes.^ 
We  have  found  that  when  X  =  (p^  —  p^)  x  trace/^/,  the 
reconstruction  algorithm  generates  the  best  results  for  PAT 
image  reconstruction. 

The  image  formation  process  described  above  is  first 
tested  using  simulated  data.  The  test  geometry  is  shown  in 
Eig.  1(a)  where  a  two-dimensional  circular  background  region 
(50.8  mm  in  diameter)  contained  four  circular  targets  (5  mm 
in  diameter  each).  The  optical  properties  for  the  background 
were  absorption  coefficient  /x^  =  0.01  mm“^  and  reduced 
scattering  coefficient  /x'  =  1.0  mm“^  while  the  optical  prop¬ 
erties  were  /x^  =  0.2  mm“^  and  /x'  =  1.0  mm“^  for  the  top 
left  and  bottom  right  targets,  and  /x^  =0.1  mm“^  and  /x' 
=  1.0  mm“^  for  the  top  right  and  bottom  left  targets.  In 
the  simulation,  a  homogeneously  distributed  area  source  was 
utilized  to  illuminate  the  whole  imaging  domain,  the  same 
as  in  our  experiments.  A  total  of  120  ultrasound  receivers 
were  equally  distributed  along  the  boundary  of  background 
region. 

Eor  the  phantom  experiments,  a  pulsed  light  from  a 
Nd:YAG  laser  (wavelength:  532  nm,  pulse  duration:  3-6  ns; 
Altos,  Bozeman,  MT)  was  sent  to  the  top  surface  of  the  cylin¬ 
drical  phantom  via  an  optical  subsystem  and  consequently 
generated  acoustic  signals.^  A  transducer  (1  MHz  central  fre¬ 
quency;  GE  Panametrics,  Waltham,  MA)  and  the  phantom 
were  immersed  in  a  water  tank.  A  rotary  stage  rotated  the 
transducer  relative  to  the  center  of  the  tank.  The  incident 
optical  fiuence  was  controlled  below  10  mJ/cm^  and  the  in¬ 
cident  laser  beam  diameter  was  5  cm.  The  complex  wave- 
field  signal  was  first  amplified  by  a  preamplifier  (gain:  17  dB, 
5  kHz  to  25  MHz;  Onda,  Corporation,  Sunnyvale,  CA),  and 
then  amplified  further  by  a  Pulser/Receiver  (GE  Panametrics, 
Waltham,  MA).  In  the  experimental  tests,  we  embedded  one 
or  two  objects  with  a  size  ranging  from  0.5  to  3  mm  in  the 
5-cm-diameter  solid  cylindrical  phantom.  The  phantom  mate¬ 
rials  used  consisted  of  Intralipid  as  scatterer  and  India  ink  as 
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Fig.  1.  Exact  distribution  of  optical  absorption  coefficient  i^a  (a)  exact  distribution  of  energy  density  (b),  the  recovered  image  using  the  new  method  (c), 
and  recovered  optical  absorption  profile  plotted  along  y  =  7.5  mm  from  the  image  shown  in  Fig.  1(c)  (d).  The  axes  (left  and  bottom)  illustrate  the  spatial  scale, 
in  mm,  whereas  the  gray  scale  (right)  records  the  normalized  absorbed  energy  density  in  arbitrary  units,  or  /x^  in  mm“^ 


absorber  with  Agar  powder  (l%-2%)  for  solidifying  the  In¬ 
tralipid  and  India  ink  solution.  We  then  immersed  the  object¬ 
bearing  solid  phantom  into  the  water  tank.  The  absorption 
of  the  background  phantom  was  0.01  mm“\  while  the  ab¬ 
sorption  coefficient  of  the  target(s)  was  0.03  mm“^.  The  re¬ 
duced  scattering  coefficients  of  the  background  medium  and 
target(s)  were  1.0  and  3.0  mm“\  respectively.  As  a  final  test, 
we  placed  human  hair-containing  phantom  (four  hairs)  into 
the  water  tank  to  show  the  high  resolution  imaging  capabil¬ 
ity  of  the  developed  scheme.  The  initial  guess  for  the  /x^  is 
0.005  mm“^  for  the  experimental  tests.  It  is  noted  that  con¬ 
stant  Gruneisen  function  T  =  c^filCp,  acoustic  velocity  (cq 
=  1495  m/s),  and  reduced  scattering  coefficient  (1.0  mm“^) 
were  employed  for  all  the  reconstructions  using  both  simu¬ 
lated  and  phantom  data. 

The  results  from  simulated  data  are  shown  in  Fig.  1,  where 
Figs.  1(a)  and  1(b)  provide  the  exact  distribution  of  /x^  and  ab¬ 
sorbed  energy  density,  respectively,  while  Fig.  1(c)  presents 
the  reconstructed  /x^  image  using  the  one-step  qPAT  algo¬ 
rithm.  We  can  see  from  Fig.  1(c)  that  the  targets  are  clearly 
identified  in  terms  of  the  target  position,  size,  and  absorption 
coefficient  value.  The  target  size  was  estimated  to  be  4.6  mm 
in  diameter  using  the  full  width  at  half-maximum  (FWHM) 
of  the  recovered  /x^  profile  in  Fig.  1(d). 


It  is  also  observed  from  Fig.  1(b)  that  the  influence  of  the 
inhomogeneous  distribution  of  optical  fiuence  on  the  energy 
density  map  is  apparent,  where  low  values  of  energy  density 
appear  in  the  center  area  of  the  four  targets  and  results  in  hol¬ 
low  phenomena  for  the  solid  targets.  However,  the  absorp¬ 
tion  coefficient  map  shown  in  Fig.  1(c)  indicates  that  the  new 
method  has  the  capability  to  completely  remove  effect  of  the 
heterogeneous  distribution  of  local  fiuence. 

The  results  from  the  first  experiment  are  shown  in  Fig.  2, 
where  Figs.  2(a)  and  2(b)  show  the  recovered  /x^  maps  us¬ 
ing  the  previous^  and  present  qPAT  methods,  respectively.  By 
estimating  the  FWHM  of  the  specific  profiles  for  the  images 
shown  in  Fig.  2(b),  we  found  the  recovered  object  sizes  were 
1.9  and  3.1  mm,  which  are  in  good  agreement  with  the  ac¬ 
tual  object  sizes  of  2.0  and  3.0  mm.  Figure  3  plots  the  recon¬ 
structed  absorption  coefficient  images  for  experimental  test  2 
using  both  the  previous^  and  present  qPAT  schemes.  The  re¬ 
covered  target  size  using  the  present  and  previous  schemes, 
as  shown  in  Fig.  3(c),  was  found  to  be  0.55  and  0.8  mm,  re¬ 
spectively,  compared  to  the  actual  target  size  of  0.5  mm.  Ob¬ 
viously  the  previous  method  overestimated  the  target  size  in 
most  cases.  We  also  found  from  Figs.  2  and  3  that  the  cur¬ 
rent  qPAT  method  is  able  to  considerably  reduce  the  effect  of 
local  optical  fiuence  and  boundary  noise.  Further,  it  is 
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Fig.  2.  Reconstructed  /Ua  image  using  the  previous  scheme  (a),  the  new  scheme  (b),  and  recovered  optical  absorption  profile  plotted  along  y  =  —7.0  mm  from 
the  images  shown  in  Figs.  2(a)  and  2(b)  (c).  The  axes  (left  and  bottom)  illustrate  the  spatial  scale,  in  mm,  whereas  the  gray  scale  (right)  records  /j.a  in  mm“'. 
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Fig.  3 .  Reconstructed  iia  images  using  (a)  the  previous,  (b)  new  scheme,  and  (c)  recovered  optical  absorption  prohle  plotted  along  =  1 .0  mm  from  the  images 
in  (a)  and  (b). 


observed  from  Fig.  4  that  based  on  the  developed  scheme, 
the  human  hairs  embedded  in  phantom  are  clearly  identi¬ 
fied  with  submillimeter  resolution  and  quantitative  absorption 
coefficient. 

In  summary,  we  have  demonstrated  experimental  evidence 
that  it  is  possible  to  directly  recover  the  absolute  /x^  image 
using  boundary  photoacoustic  measurements  coupled  with 
the  photon  diffusion  equation  in  just  one  step.  The  method 
described  is  able  to  quantitatively  reconstruct  absorbing  ob¬ 
jects  with  different  sizes  and  optical  contrast  levels.  In  par¬ 
ticular,  it  is  noted  that  this  method  does  not  depend  on  any 
calibration  procedure  because  it  allows  for  the  use  of  rela¬ 
tive  incident  laser  source  strength  and  normalized  boundary 
measurements  of  acoustic  pressure.  The  source  intensity  and 
boundary  condition  coefficient  can  be  optimized  and  localized 
based  on  the  photon  diffuse  equation  by  minimizing  /  with 
the  normalized  acoustic  measurements  p:  f  =  YlfLi 
—  pi)^.  As  such,  the  reconstruction  of  absorption  coefficient 
with  the  new  algorithm  does  not  depend  on  the  absolute  val¬ 


ues  of  absorbed  energy  density  and  optical  fiuence  whereas 
the  previous  methods  are  all  dependent  on  these  values.  Ad¬ 
ditionally,  we  note  that  different  transducer  sensitivity  may 
lead  to  a  different  scaling  factor  of  the  reconstructed  absorp¬ 
tion  coefficient  images.  Routine  calibrations  on  the  transducer 
response  had  to  be  done  in  our  previous  study. ^  However,  in 
the  new  scheme  presented  here,  we  use  normalized  acoustic 
data  in  arbitrary  units  (i.e.,  the  largest  value  of  pressure  is  1.0) 
to  recover  the  absorption  coefficient,  reducing  the  impact  of 
transducer  sensitivity  significantly.  It  should  be  pointed  out, 
however,  that  the  choice  of  the  initial  guess  of  /x^  does  have 
some  effect  on  the  recovered  results.  Finally,  further  investi¬ 
gation  is  warranted  to  consider  optical  scattering,  inhomoge¬ 
neous  acoustic  properties,  inhomogeneous  light  illumination 
field,  and  multispectral  photoacoustic  measurements  within 
the  framework  of  this  new  method. 
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While  brain  imaging  and  electrophysiology  play  a  central  role  in  neuroscience  research  and  in  the  evaluation 
of  neurological  disorders,  a  single  noninvasive  modality  that  offers  both  high  spatial  and  temporal  resolution 
is  currently  not  available.  Here  we  show  in  an  acute  epilepsy  rat  model  that  photoacoustic  tomography  (PAT) 
can  noninvasively  track  seizure  brain  dynamics  with  both  high  spatial  and  temporal  resolution,  and  at  a 
depth  that  is  clinically  relevant.  The  noninvasive  yet  whole  surface  and  depth  capabilities  of  the  PAT  system 
allowed  us  to  actually  see  what  is  happening  during  ictogenesis  in  terms  of  seizure  onset  and  spread.  Both 
seizure  onset  and  propagation  were  tomographically  detected  at  a  spatial  resolution  of  150  pm  and  a  tempo¬ 
ral  resolution  of  300  ms,  respectively.  The  current  study  lends  support  to  the  theory  that  seizure  onset  and 
spread  involves  a  rich  interplay  between  multiple  cortical  and  subcortical  brain  areas  during  the  onset  and 
spread  of  epileptic  seizures.  Dynamical  changes  of  vasculature  during  epileptiform  events  were  also  detected 
with  high  spatiotemporal  resolution.  Together,  these  findings  suggest  that  PAT  represents  a  powerful  tool  for 
noninvasively  mapping  seizure  onset  and  propagation  patterns,  and  the  ‘functional’  connectivity  within 
epileptic  brain  networks. 

©  2012  Elsevier  Inc.  All  rights  reserved. 


Introduction 

Epilepsy  is  a  common,  chronic,  neurological  disorder  character¬ 
ized  by  seizures.  Three  percent  of  people  will  be  diagnosed  with 
epilepsy  at  some  time  in  their  lives  (Hauser  et  al.,  1996).  Indeed,  ap¬ 
proximately  50  million  people  worldwide  have  epilepsy,  and  20  to 
30%  of  these  patients  are  refractory  to  all  forms  of  medical  treatment 
(Hauser,  1992).  In  most  cases,  seizures  are  controlled,  although  not 
cured,  with  anticonvulsant  medication.  Seizure  types  are  classified 
according  to  whether  the  source  of  the  seizure  is  localized  (partial 
or  focal  onset  seizures)  or  distributed  (generalized  seizures) 
(Greenfield  et  al,  2011).  Localization-related  epilepsies  arise 
from  an  epileptic  focus.  A  partial  seizure  may  spread  within  the 
brain,  a  process  known  as  secondary  generalization.  Generalized 
epilepsies,  in  contrast,  arise  from  many  independent  foci  (multifocal 
epilepsies)  or  from  epileptic  circuits  that  involve  the  whole  brain.  In 
epilepsies  of  unknown  localization,  it  remains  unclear  as  to  whether 
they  arise  from  a  portion  of  the  brain  or  from  more  widespread 
circuits. 
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For  those  patients  with  medically  intractable  focal  epilepsy,  the 
best  treatment  option  is  resective  brain  surgery  (Berg  et  al.,  2007; 
Birbeck  et  al.,  2002;  Duncan  et  al.,  2006).  Although  several  factors 
can  impact  the  success  of  epilepsy  surgery,  the  primary  reason  of  fail¬ 
ure  is  the  incomplete  mapping  of  the  local  epilepsy  network  which 
results  in  incomplete  resection  of  epileptogenic  foci  (Engel,  2004; 
Jeha  et  al.,  2007).  Much  of  the  attention  on  epilepsy  surgery  has 
been  directed  at  identifying  single  neuronal  populations.  This  ap¬ 
proach  has,  in  many  cases,  led  to  failed  surgical  outcomes,  because 
seizures  typically  involve  groups  of  neurons  interacting  both  locally 
and  across  several  cortical  and  subcortical  brain  regions.  A  better  un¬ 
derstanding  of  the  regional  interactions  occurring  at  the  site  of  sei¬ 
zure  onset  and  spread  may  provide  important  insights  about  the 
pathophysiology  of  seizures  and  aid  with  accurate  brain  mapping 
and  resection  of  the  epileptic  focus. 

In  theory,  removing  the  focus  should  result  in  a  patient's  seizures 
being  cured.  However,  there  is  much  evidence  to  suggest  that  the 
focus  is  more  of  a  region  of  seizure  onset  with  a  number  of  sites 
that  can  act  independently  to  initiate  seizures  (Thom  et  al,  2010a, 
b).  Seizures  in  animal  models  and  in  people  often  have  a  multifocal 
or  broadly  synchronized  onset.  The  best  evidence  for  multifocality 
within  the  seizure  onset  zone  comes  from  surgical  experience  with 
intracranial  monitoring.  The  challenge  of  mapping  the  epileptic 
focus  stems  from  the  observation  that  the  pathology  associated  with 
focal  epilepsy  is  often  distributed  across  a  number  of  brain  sites 
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(Bertram,  2009)  and  that  current  diagnostics  methods  frequently  fall 
short  of  identifying  such  sites.  Animal  studies  indicate  that  the  neu¬ 
rons  involved  in  the  epileptic  circuitry  have  enhanced  excitability 
throughout  (Bertram  et  al,  1998;  Fountain  et  al,  1998;  Mangan  et 
al,  2000).  The  implication  of  these  observations  is  that  each  of 
the  sites  could  act  independently  to  initiate  a  seizure  or,  potential¬ 
ly,  to  drive  another  site  into  a  seizure.  Thus,  in  focal  epilepsy,  one 
may  view  a  cortical  region  as  a  broad  seizure  onset  zone,  with 
the  potential  that  multiple  foci  can  act  as  a  seizure  focus  for  any 
given  seizure. 

Much  of  our  understanding  about  focal  seizure  circuitry  comes  from 
electrophysiological  recording  methods.  Although  electrophysiology  is 
currently  the  ‘gold  standard’  in  mapping  the  epileptic  focus,  it  is  often  in¬ 
adequate  to  define  the  boundary  of  the  epilepsy  circuitry  due  to  spatial 
sampling  limitations  and  volume  conduction.  What  we  truly  desire  is  a 
brain  mapping  modality,  which  could  give  a  high-resolution  real  time 
spatial  and  temporal  ‘read  out’  of  the  dynamics  of  cortical  processing 
and  seizures.  It  has  been  well  established  that  optical  contrast  is  highly 
sensitive  to  neuronal  activity  (Grinvald  et  al.,  1988;  Hill  and  Keynes, 
1949).  The  high  optical  contrast  is  largely  due  to  the  changes  in  blood 
volume  and  blood  oxygenation,  both  of  which  substantially  increase 
during  seizures,  and  increase  in  oxygenation  metabolic  rate  results  in  in¬ 
creased  demands  on  autoregulatory  mechanisms  (Folbergrova  et  al., 
1981 ).  Taking  the  advantage  of  the  oxygenation  dependence  in  the  op¬ 
tical  absorption  spectrum  of  hemoglobin,  optical  imaging  or  intrinsic  op¬ 
tical  signal  (lOS)  methodologies  have  provided  excellent  surface  maps 
of  epileptic  focus  (Haglund  and  Hochman,  2004;  Haglund  et  al.,  1992). 
Although  near-infrared  spectroscopy  (NIRS)  (Steinhoff  et  al.,  1996, 
Sokol  et  al.,  2000,  Watanabe  et  al.,  2000)  and  blood-oxygen-level-de- 
pendent  (BOLD)  MRI  can  assess  oxygen  saturation  with  endogenous 
contrasts,  BOLD  MRI  is  sensitive  only  to  HbR  with  low  temporal  resolu¬ 
tion  (Ogawa  et  al.,  1990).  In  addition,  the  major  limitation  of  NIRS  and 
lOS  is  that  each  provides  only  surface  depth  information.  Greater 
depth  information  can  be  obtained  by  using  tomographic  reconstruc¬ 
tion  methods  such  as  diffuse  optical  tomography  (Bluestone  et  al., 
2001 ;  Boas  et  al.,  2001 ;  Jiang  et  al.,  1 996;  Oleary  et  al.,  1 995).  Whereas 
diffuse  optical  tomography  has  low  spatial  resolution,  laser-induced 
photoacoustic  tomography  (PAT)  has  both  superior  spatial  and  tempo¬ 
ral  resolution.  PAT  detects  absorbed  photons  ultrasonically  by 
employing  the  photoacoustic  effect.  It  combines  both  high  contrast 
and  spectroscopic  specificity  based  on  the  optical  absorption  of  both 
oxy-  and  deoxy-hemoglobin  with  high  ultrasonic  spatial  resolution 
(Kruger  and  Liu,  1994;  Laufer  et  al.,  2007;  Xiang  et  al.,  2007;  Yuan  et 
al.,  2007).  Relative  to  other  optical  imaging  modalities,  PAT  has  the  ad¬ 
vantage  of  mitigating  both  scalp  and  skull  light  scattering  by  a  factor  of 
-1000.  The  end  result  is  that  PAT  allows  for  high  spatial  resolution 
imaging  of  brain  at  a  depth  considerably  beyond  the  soft  depth  limit 
of  conventional  optical  imaging  techniques  such  as  confocal  microscopy 
(Sipkins  et  al.,  2005),  two-photon  microscopy  (Denk  et  al.,  1990),  and 
optical  coherence  tomography  (Huang  et  al.,  1991 ).  The  strong  prefer¬ 
ential  optical  absorption  of  hemoglobin  makes  PA  imaging  to  have  a 
better  imaging  contrast  than  ultrasound  (US) ;  as  it  can  be  difficult  to  vi¬ 
sualize  the  microvessels  with  pulse-echo  US  owing  to  the  weak 
echogenicity  (Mace  et  al.,  2011 ). 

In  this  study,  PAT  was  employed  to  image  seizures  in  an  experimen¬ 
tal  acute  bicuculline  methiodide  model  of  focal  epilepsy.  Bicuculline  is  a 
light-sensitive  competitive  antagonist  of  GABAa  receptors  that  mimics 
focal  epilepsy  when  applied  to  brain  tissue.  During  focal  application  of 
bicuculline  into  the  brain  cortex,  brains  were  imaged  noninvasively 
with  a  novel  PAT  system  that  has  three  orders  of  magnitude  higher  tem¬ 
poral  resolution  and  four-fold  higher  spatial  resolution  relative  to  our 
previous  PAT  prototype  (Zhang  et  al.,  2008).  Off-line,  we  employed 
measures  of  brain  connectivity  to  further  identify  the  functional  anato¬ 
my  implicated  in  focal  cortical  seizures.  The  high  spatial  and  temporal 
sampling  of  the  novel  PAT  system  allowed  for  the  first  time  the  com¬ 
plete  mapping  of  an  epileptiform  event  in  vivo. 


Methods 

Animals 

Male  Sprague-Dawley  rats  (Harlan  Labs,  Indianapolis,  IN)  weighing 
50-60  g  on  arrival  were  allowed  one  week  to  acclimate  to  the  12-h 
light/dark  cycle  and  given  food  and  water  ad  libitum.  All  procedures 
were  approved  by  the  University  of  Florida  Animal  Care  and  Use  Com¬ 
mittee  and  conducted  in  accordance  with  the  National  Institutes  of 
Health  Guide  for  the  Care  and  Use  of  Experimental  Animals. 


PAT  imaging 

Light  from  a  Ti:  Sapphire  laser  tunable  (690  to  950  nm)  was  delivered 
through  the  skull  to  the  brain  through  an  optical  fiber  (Fig.  la).  The 
energy  of  each  laser  pulse  was  detected  by  a  photodiode  for  calibration. 
A  192  element  full-ring  transducer  array  was  used  to  capture  the 
photoacoustic  (PA)  signals  generated  by  the  laser  light.  The  192  channel 
data  acquisition  system  consisted  of  preamplifiers,  secondary  stage  am¬ 
plifiers  (for  optimizing  the  signal-to-noise  ratio),  and  a  3:1  electronic 
multiplexer  coupled  with  a  64-channel  analog-to-digital  converter. 
Each  ultrasonic  detector  had  a  5-MHz  central  frequency  and  a  70% 
nominal  bandwidth  with  a  diameter  of  6  mm  (Blatek,  Inc,  PA,  USA). 


Fig.  1.  Real  time  PAT  system  for  seizure  dynamics,  (a.)  Schematic  of  the  real  time  PAT 
system.  A  192-element  full-ring  transducer  array  was  used  to  capture  the  PA  signal 
during  seizure  onset,  (b.)  Noninvasive  PAT  imaging  of  a  rat  brain  in  vivo  with  the 
skin  and  skull  intact.  MCA,  middle  cerebral  artery:  RH,  right  hemispheres:  LH,  left 
hemispheres:  LOB,  left  olfactory  bulbs:  ROB,  Right  olfactory  bulbs,  (c.)  Open-skull  pho¬ 
tograph  of  the  rat  brain  surface  acquired  after  the  PAT  experiment.  Numbers  1-5  indi¬ 
cate  the  corresponding  blood  vessels  in  the  PAT  image  and  rat  brain  photograph. 
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Electrode  implantation  surgery 

Animals  were  anesthetized  intraperitoneally  with  1  g/kg  of  body 
weight  dose  of  urethane.  Two  300  |um  diameter  stainless  steel  screw 
electrodes  were  implanted  within  the  skull  for  obtaining  subdural 
multichannel  cortical  local  field  potentials  data  ( —  0.3  mm  posterior, 
3  mm  lateral  (right)  of  bregma,  1  mm  ventral)  based  on  coordinates 
from  a  rat  brain  atlas  (Paxinos,  1998).  One,  300  |um  diameter  stainless 
steel  screw  electrode  (FHC,  Bowdoin,  ME)  was  implanted  as  a  refer¬ 
ence  electrode  into  the  midline  occipital  bone.  Cortical  local  field  po¬ 
tentials  were  obtained  using  a  Tucker  Davis  Pentusa  (Tucker  Davis 
Technologies,  Alachua,  FL)  neural  recording  system  at  12  kHz,  digi¬ 
tized  with  16  bits  of  resolution,  and  band  pass  filtered  from  0.5  to 
6  kHz. 

Induction  of  seizures 

Rats  (n  =  10)  received  10  pi  of  1.9  mM  bicuculline  methiodide 
and  10  pi  normal  saline  into  the  left  and  right  parietal  cortex,  re¬ 
spectively.  The  infusion  was  performed  through  the  previously 
implanted  electrode  sites  at  a  rate  of  0.3  pl/min.  The  infusion  sys¬ 
tem  consisted  of  a  100  pi  gas-tight  syringe  (Hamilton,  Reno,  NV) 
driven  by  a  syringe  pump  (Cole-Parmer,  Vernon  Hills,  IL)  connected 
to  polyaryletheretherketone  (PEEK)  tubing  (ID  =  0.381  mm,  OD  = 
0.794  mm,  length -0.5  m,  Upchurch  Scientific,  Oak  Harbor,  WA). 
The  PEEK  tubing  was  coupled  to  a  silica  cannula  (ID  =  50  pm, 
OD  =  147  pm.  Polymicro  Technologies,  Phoenix,  AZ)  via  a  microfiuidic 
connector.  Cortical  local  field  potentials  were  recorded  5  min  before 
each  injection  and  continued  for  up  to  30  min  thereafter. 

Image  analysis 

A  custom  software  utility  was  written  and  incorporated  into  the 
computer  software  package  MATLAB  (MathWorks,  Massachusetts, 
USA)  to  analyze  and  display  recorded  data.  This  software  enabled 
the  reconstruction  of  the  PAT  images  and  the  determination  of  the 
blood  vessel  diameter.  Amira  (version  5.3.3,  TGS  Template  Graphics 
Software)  was  used  for  three-dimensional  reconstruction  of  the  sei¬ 
zure  foci.  We  used  the  frequency  domain  Granger  causality  method¬ 
ology  to  evaluate  the  dynamic  interactions  within  the  seizure 
circuitry  and  the  associated  brain  networks  (Brovelli  et  al.,  2004; 
Granger,  1969;  Wang  et  al.,  2007).  The  time-series  data  were  selected 
from  200  sets  of  PAT  images  sampled  at  3  Hz  within  the  seizure  onset 
period,  and  at  each  time  point,  we  chose  9  (regions  of  interest)  ROI 
(R1-R9)  and  then  the  PA  signal  was  averaged  within  each  ROI,  as 
shown  in  Fig.  5bl.  Changes  in  photoacoustic  signals  were  quantified 
as  —  AA/A  in  Fig.  5b3,  right.  The  PA  data  shown  in  Figs.  3-5  were 
not  averaged.  The  image  color  scale  was  determined  by  the 
photoacoustic  signal  intensity  in  arbitrary  units  (between  0  and  1). 

Results 

Noninvasive  epileptic  foci  localization 

Temporal  and  spatial  resolutions  of  the  PAT  system  were  experi¬ 
mentally  determined.  PA  images  depicting  frames  from  a  10-second 
sequence  of  dynamic  ink  fiow  through  a  0.3  mm  diameter  tube  was 
shown  in  Movie  SI.  The  movie  shows  that  the  temporal  resolution 
of  this  imaging  system  was  -0.33  s/frame.  Cross-sectional  PA  image 
of  a  two-copper-rods  phantom  was  used  to  test  the  spatial  resolution 
of  the  PAT  imaging  system.  The  two  copper  rods  (diameter:  0.05  mm) 
were  embedded  in  a  cylindrical  tissue-mimicking  phantom  at  a  depth 
of  10  mm.  The  center-to-center  distance  between  the  two  copper 
rods  was  approximately  2.2  mm.  The  spatial  resolution  of  the  imag¬ 
ing  system  calculated  with  Rayleigh's  law  was  1 50  pm.  The  spatial 


resolution  of  the  PAT  system  was  inversely  related  to  the  bandwidth 
of  the  ultrasonic  transducer. 

The  noninvasive  PAT  image  of  the  rat  cortical  vasculature  (Fig.  lb) 
matched  well  with  an  anatomical  atlas  (Paxinos,  1998)  (Fig.  Ic) 
obtained  after  the  PA  imaging.  The  brain  structures  including  the 
middle  cerebral  artery,  right  hemispheres,  left  hemispheres,  left  ol¬ 
factory  bulbs,  and  right  olfactory  bulbs  are  clearly  shown  in  the  PAT 
image.  Numbers  1-5  marked  in  the  image  indicate  that  the 
micro-blood  vessels  with  a  diameter  of  less  than  100  pm  are  also 
seen,  which  again  correspond  well  with  the  rat  brain  photograph. 
This  allowed  for  imaging  objects  of  50  pm  in  diameter  with  a  spatial 
resolution  of  150  pm  (Figs.  2b,  c). 

Thirty  minutes  following  the  microelectrode  placement,  PAT  imag¬ 
ing  was  performed  to  generate  baseline  tissue  absorption  maps,  visual¬ 
ize  micro  blood  vessels,  and  morphological  cortical  landmarks  (Figs  3b, 
c).  Subsequently,  rats  (n  =  10)  received  10  pi  of  1.9  mM  bicuculline 
methiodide  and  10  pi  normal  saline  into  the  left  [— l.OAP,  2.5ML, 
2.4DV]  and  right  [  —  l.OAP,  2.5  ML,  0.4DV]  parietal  cortex,  respectively, 
based  on  coordinates  from  a  rat  brain  atlas  (Paxinos,  1998).  Immediate¬ 
ly  following  the  focal  infusions,  PAT  imaging  was  repeated  to  visualize 
changes  in  the  injection  site  tissue  absorption.  An  increase  in  tissue  ab¬ 
sorption  was  observed  in  the  bicuculline  methiodide  cortical  injection 
site  and  surrounding  region  (Fig.  3c,  left),  but  not  in  the  saline  injection 
site  or  surrounding  region  (Fig.  3c,  right).  PAT  scanning  for  subcortical 
changes  was  also  performed.  Slices  5,  7,  and  9  (Fig.  3d)  are  the  images 
obtained  3,  5  and  7  mm  below  the  scalp,  respectively.  Fig.  3e  is  the 
three  dimensional  rendering  of  the  epileptic  foci  obtained  from  differ¬ 
ent  tomographic  layers  (movie  S2).  Each  image  exhibited  the  variable 
patterns  of  seizure  onset  and  propagation  both  at  the  cortical  and  sub¬ 
cortical  regions.  Seizures  were  confirmed  by  concomitant  time-locked 
PAT/video-electroencephalography  (Fig.  3a).  Experiments  were  repeat¬ 
ed  for  each  animal  at  2-hour  intervals. 

Real  time  monitoring  of  epileptic  events 

Epileptiform  events  were  recorded  with  PAT  from  a  focal  region  of 
interest  of -2x3  mm  (Fig.  4a).  We  observed  a  significant  optical  ab¬ 
sorption  change  directly  associated  with  de-oxygenation  at  a  wave¬ 
length  of  755  nm  (Fig.  4c).  Furthermore,  at  the  seizure  onset  at  1  min 
6.267  s,  the  seizure  focus  measured  0.2053  mm^  and  increased  to 
0.5004  mm^  as  the  electrographic  seizure  time-series  increased  in  fre¬ 
quency  and  amplitude  (Fig.  4b).  Corresponding  rate  changes  (0.33  s) 
spike  and  wave  discharges  and  PA  images  were  observed  in  each  exper¬ 
imental  trial.  The  seizure  onset  dynamics  captured  photoacoustically 
are  best  appreciated  in  movies  displaying  how  the  seizure  was  generat¬ 
ed  and  varied  over  time  (Movie  S3).  PAT  images  suggest  that  in  addition 
to  ictal  spread  from  the  primary  focus,  homotopic  foci  are  seen  in  the 
contralateral  parietal  cortex  (Fig.  5al-4).  The  PA  signal  in  the  contralat¬ 
eral  foci  was  smaller  in  magnitude  and  delayed  in  time  compared  to  the 
signal  recorded  from  the  primary  seizure  foci.  We  also  observed  a  re¬ 
gion  of  inverted  optical  signal  immediately  surrounding  the  seizure 
onset  zone  (Fig.  5b2-3).  We  performed  analysis  of  the  optical  signal  at 
a  distance  of  at  least  2  mm  (±0.01)  away  from  the  edge  of  the  focus. 
Our  results  demonstrate  that  the  photoacoustic  signal  was  inversely  re¬ 
lated  to  the  optical  signal  recorded  from  the  focus  (Fig.  5b3). 

To  evaluate  the  dynamic  interactions  within  the  seizure  circuitry 
and  capture  the  associated  networks  we  used  the  frequency  domain 
Granger  causality  (Brovelli  et  al.,  2004;  Granger,  1969;  Wang  et  al., 
2007).  PAT  time  series  data  was  selected  from  200  sets  of  PAT  images 
sampled  at  3  Hz  within  the  seizure  onset  period.  At  each  time  point 
we  chose  9  ROIs  (Fig.  5bl)  and  used  the  averaged  optical  absorption 
within  each  region  for  the  network  analysis.  Nine  sets  of  time  series 
analysis  were  classified  into  3  groups  including  the  primary  ictal  onset 
area  and  corresponding  4  surrounding  regions  of  interest  including 
the  primary  or  initiating  seizure  focus,  contralateral  homotopic  foci, 
and  2  regions  of  interest  surrounding  the  primary  focus  in  the 
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Fig.  2.  Temporal  and  spatial  resolution  demonstration  of  the  PAT  system,  (a.)  Temporal  resolution  demonstration  of  the  PAT  system.  Photoacoustic  images  depicting  frames  from  a  10-second 
sequence  of  dynamic  ink  flow  through  a  0.3  mm  diameter  tube  (Movie  SI ).  The  movie  shows  that  the  temporal  resolution  of  this  imaging  system  was  -0.33  s/frame.  (b.)  Cross-sectional  PA 
image  of  a  two-copper-rods  phantom.  The  two  copper  rods  (diameter:  0.05  mm)  were  embedded  in  a  cylindrical  tissue-mimicldng  phantom  at  a  depth  of  10  mm.  The  center-to-center  dis¬ 
tance  between  the  two  copper  rods  was  approximately  2.2  mm.  (c.)  The  normalized  profile  of  the  reconstructed  PA  image  shown  in  (b.)  along  y=  10  mm.  The  half-  and  quarter-amplitude 
widths  were  obtained  by  measuring  the  distance  between  points  AB  and  between  CD,  respectively.  The  spatial  resolution  of  the  imaging  system  calculated  with  Rayleigh's  law  was  150  pm. 


contralateral  cortex.  The  mi(idle  cerebral  artery,  primary  seizure  focus, 
and  secondary  homotopic  seizure  focus  were  also  analyzed.  Results 
from  group  1  analysis  demonstrate  the  causal  influence  from  the  prima¬ 
ry  focus  to  the  surrounding  regions  of  interest  (Fig.  5bl ).  The  negative 
zero-lag  correlation  (  —  0.68)  between  both  the  primary  and  secondary 
foci  showed  that  neural  activities  within  these  two  regions  changed  in 
opposite  direction  over  time.  Fig.  5b4  shows  the  Granger  causality 
analysis  for  group  2  data  where  we  see  that  the  primary  and  secondary 
contralateral  brain  foci  influenced  each  other  mutually  with  a  stronger 
influence  of  the  primary  focus  relative  to  the  contralateral  focus.  More¬ 
over,  we  suggest  that  the  primary  seizure  focus  may  influence  the  con¬ 
tralateral  foci  as  previously  suggested  using  optical  intrinsic  imaging  in 
epilepsy  (Khalilov  et  al.,  2003)  and  microelectrode  recordings  of  hippo¬ 
campus  local  field  potentials  (Cadotte  et  al.,  2010).  Finally,  from  group  3 
analyses  (Fig.  5b5)  we  found  that  the  primary  focus  strongly  influenced 
the  hemodynamic  change  in  the  middle  cerebral  artery,  which  in  turn 
showed  influence  on  the  contralateral  homotopic  foci. 

Dynamical  changes  of  vasculature  during  interictal  discharges 

Fig.  6a  shows  image  of  the  rat  cortex  vasculature  through  the  in¬ 
tact  scalp  and  skull.  The  red  arrow  indicates  the  microvascularity 
along  the  direction  marked  by  the  yellow  dashed  line,  representing 
a  cross-sectional  scanning  using  a  50  MHz  ultrasound  transducer. 
Positive  and  negative  acoustic  peaks  induced  by  a  25  jum  blood  vessel 


were  sorted  out  as  target  signals  to  track  the  change  of  the  vessel  diam¬ 
eter.  Typical  PA  signals  from  the  targeted  vessel  at  different  times  show 
apparent  vessel  vasomotion  by  -40%  (mean  =  40%,  P<0.03,  N  =  20) 
during  the  interictal  discharges  (Fig.  6b).  Local  field  potential  recordings 
(Fig.  6c)  showed  that  interictal  spikes  had  a  strong  correlation  with  the 
discrete  change  in  vessel  size  observed  photoacoustically  (Fig.  6d;  error 
bars  ±  s.d.  was  calculated  from  20  consecutive  spikes). 

Discussion 

The  main  finding  is  that  PAT  imaging  is  a  novel  tool  for  noninva- 
sively  mapping  seizure  dynamics  with  both  high  spatial  and  temporal 
resolution.  This  study  investigated  the  hemodynamic  changes  during 
focal  seizure  onset  and  propagation  in  an  acute  model  of  focal  epilep¬ 
sy.  The  experimental  results  suggest  that  the  increase  in  local  and 
surrounding  brain  tissue  absorption  was  due  to  the  focal  bicuculline 
methiodide  induced  seizure  rather  than  other  factors  such  as  stab 
cortex  tissue  wound  injection.  Epileptiform  events,  including 
interictal  spikes,  ictal  onset  and  spread  were  identified  in  near 
real  time.  To  the  best  of  our  knowledge,  this  is  the  first  experimen¬ 
tal  evidence  for  millisecond  temporal  resolution,  centimeters  depth, 
and  micron  scale  imaging  of  focal  seizure  onset  and  spread  by  a 
noninvasive  technique. 

Quantitative  photoacoustic  image  reconstruction  can  be  used  to 
resolve  the  light  scattering  problem  with  the  deeper  penetration  in 
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Fig.  3.  Noninvasive  epileptic  foci  localization,  (a.)  EEG  recordings  of  seizure  onset  after  BMl  injection,  (b.)  Schematic  showing  the  location  of  BMl  injection  and  saline  injection  (control).  B,  BMl 
injection:  S,  saline  injection,  (c.)  Reconstructed  PAT  image  before  (left)  and  after  (right)  the  BMl  injection.  The  BMl  and  saline  were  injected  into  the  left  and  right  parietal  neocortex,  respec¬ 
tively.  A  significant  increase  of  optical  absorption  is  seen  in  the  region  of  the  BMl  injection,  while  no  absorption  contrast  is  observed  in  the  region  of  the  saline  injection  relative  to  its  surround¬ 
ings.  (d.)  Series  of  PAT  images  from  three  representative  transver  planes  parallel  to  the  skin  surface  with  755  nm  wavelength.  Slices  5, 7,  and  9  are  the  images  obtained  3, 5  and  7  mm  under  the 
skin,  respectively.  The  images  show  different  spatial  patterns  at  the  seizure  foci  at  different  tomographic  layers,  (e.)  Three  dimensional  rendering  of  the  epileptic  foci  from  different  tomograph¬ 
ic  layers  (Movie  S2). 


the  tissue  based  on  compensating  the  reconstructed  image  with  a 
pre-calculated  optical  fluence  distribution  (Yao  et  al,  2009;  Yuan  et 
al,  2007).  It  is  possible  to  achieve  higher  resolution  at  cellular  and 
even  sub-cellular  levels  to  image  the  neuronal  circuitry  in  vivo  by 
use  of  photoacoustic  microscopy  at  a  depth  of  more  than  3  mm 
(Shelton  and  Applegate;  Hu  et  al.,  2010).  The  findings  imply  that 
PAT  has  potential  for  noninvasive  real  time  brain  mapping  of  cortical 
processing  and  seizures. 

A  key  advantage  of  our  noninvasive  PAT  system  for  epileptiform 
event  monitoring  is  the  real  time  imaging  ability.  Seizures  are  gener¬ 
ated  by  complex  interactions  of  a  large  group  of  neurons  in  which  the 
population  of  neurons  involved  varies  largely  from  moment  to  mo¬ 
ment  (Schwartz  and  Bonhoeffer,  2001).  By  achieving  1000  times 


faster  data  acquisition  of  the  current  PAT  imaging  system  relative  to 
our  previous  PAT  system  (Zhang  et  al.,  2008),  we  were  able  to  observe 
that  the  spatial  aspects  of  the  seizure  focus  varied  over  time  as 
suggested  by  the  changes  in  hemodynamics  (Fig.  4c).  At  present, 
the  imaging  speed  was  limited  only  by  the  10  Hz  laser  pulse- 
repetition-rate.  Since  considerably  faster  lasers  of  -500  Hz  pulse- 
repetition-rates  are  now  commercially  available,  the  temporal  resolu¬ 
tion  of  PAT  will  soon  reach  a  level  of  a  few  milliseconds.  The  light 
fluence  on  the  skin  was  less  than  5  mj/cm^,  well  within  the  American 
National  Standards  Institute  (Institute,  2007)  safety  limits.  This, 
coupled  with  the  compactness  of  an  optical  system  makes  it  possible 
to  bring  a  PAT  system  to  the  bedside  for  real  time  chronic  noninvasive 
monitoring  of  seizures  (Yang  et  al,  2007). 
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Fig.  4.  Real-time  monitoring  of  ictai  onset,  (a.)  A  PA  image  showing  the  cerebral  vasculature  and  the  position  of  BMl  injection.  ROl,  region  of  interest;  Scale  bar:  2  mm.  (b.)  EEG 
recordings  showing  the  seizure  onset  about  1  min  after  BMl  injection,  (c.)  PAT  images  recording  the  ictai  onset  in  real  time.  At  1  min  6.267  s  the  area  of  the  focus  was 
0.2053  mm^.  The  area  of  the  focus  then  increased  in  size  over  the  next  2  s  to  0.5004  mm^,  corresponding  to  an  increase  in  the  amplitude  of  the  EEG  spikes.  The  area  of  the  seizure 
focus  was  derived  from  the  PAT  images  by  thresholding  to  a  pixel  value  one  standard  deviation  above  the  pixel  values  from  the  area  of  the  focus  during  the  control  conditions.  Scale 
bar:  500  pm  (Movie  S3). 


One  hypothesis  is  that  long-standing  epilepsy  may  lead  to  the  de¬ 
velopment  of  secondary  epileptogenic  regions  located  at  a  site  distant 
from  the  original  focus,  a  factor  that  may  reduce  the  likelihood  of  suc¬ 
cessful  epilepsy  pre-surgical  planning  (Cendes  et  al,  1995).  The  abil¬ 
ity  to  identify  the  epileptic  focus  and  associated  network  may  lead  to 
better  epilepsy  surgery  outcomes  for  many  individuals.  Various 
methods  such  as  cross  correlation  and  Granger  causality  have  been 
explored  for  assessing  the  dynamic  directional  relationships  among 
brain  regions  using  time  series  data  (Bressler  et  al.,  2008;  Kaminski 
et  al.,  2001).  We  used  Granger  causality  as  a  way  to  quantify  the 
dynamic  interactions  amongst  the  primary  seizure  focus,  secondary 
seizure  foci,  and  middle  cerebral  artery  based  on  the  time-series 
photoacoustic  data.  Specifically,  Granger  causality  allowed  for  assess¬ 
ment  of  the  magnitude  and  direction  of  temporal  relationships  during 


overlapping  PAT  time-series  windows.  Results  further  exemplified 
the  temporal  aspects  of  seizure  onset,  seizure  secondary  spread,  and 
seizure  termination.  Collectively,  PAT  image  analysis  and  tools  for  ef¬ 
fective  connectivity  may  result  in  a  better  understanding  of  epileptic 
networks  at  high  spatiotemporal  resolution. 

Little  data  are  available  regarding  a  link  between  vascular  changes 
and  seizures  (Ndode-Ekane  et  al.,  2010).  We  found  that  the  vasomo¬ 
tor  phenomena  of  micro  blood  vessels  correlate  with  interictal  spike 
and  wave  discharges  (Fig.  6d).  Interictal  spikes  generated  a  strong  ce¬ 
rebral  metabolic  change  that  induced  cerebral  vessel  dilation  and  a 
focal  increase  in  cerebral  blood  flow  provided  oxygenated  hemoglo¬ 
bin  to  hypermetabolic  neurons,  with  a  corresponding  increase  in 
total  hemoglobin  (Ma  et  al,  2009).  Vasomotion  during  seizure  onset 
and  Spread  was  clearly  caused  by  the  oscillation  in  vessel  diameter 
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Fig.  5.  In  vivo  maps  of  mirror  foci  propagation,  (a.)  Left,  the  PA  image  of  the  rat  brain;  the  yellow  dotted  rectangular  shows  the  ROl  for  analysis;  Right,  noninvasive  imaging  of  the  primary 
and  mirror  foci  in  bilateral  homotopic  cortex,  al-4  shows  how  mirror  foci  were  generated  and  diminished.  White  arrow  shows  the  location  of  mirror  foci.  The  PA  signal  in  the  mirror  foci 
was  smaller  in  magnitude  and  delayed  in  time  compared  to  the  signal  recorded  from  the  primary  foci,  (b.)  Granger  analysis  of  seizure  foci,  mirror  foci  and  middle  cerebral  artery;  bl,  a  PAT 
image  showing  the  nine  selected  ROls;  b2.  Granger  analysis  of  the  image  around  the  foci  and  its  ‘surround’  area;  b3,  PA  signal  from  the  foci  (red)  and  surround  (blue)  during  an  ictal  event. 
The  PA  signal  surrounding  the  foci  (blue  line)  shows  inverted  signal  compared  to  that  from  the  seizure  foci  (red  line)  corresponding  to  neuronal  inhibition.  b4,  connectivity  relation  be¬ 
tween  the  primary  and  mirror  focus;  b5,  connectivity  relation  between  the  primary  focus,  mirror  focus  and  middle  cerebral  artery.  The  line  width  stands  for  the  relative  Granger  causal 
influence  strength.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  of  this  article.) 


and  the  expansion  of  the  vessel  cross  section.  We  reasoned  that  the 
dilation  resulted  in  an  active  cerebral  microvasculature  response  to 
increased  metabolic  activity  accompanying  the  seizure  (Myers  and 
Intaglietta,  1976).  The  combination  of  multi-scale  imaging  ability  and 
endogenous  hemoglobin  contrast  makes  PAT  a  promising  tool  for  imag¬ 
ing  microvasculature  during  seizure  onset.  Future  quantitative  PAT 
studies  will  allow  us  to  resolve  many  of  the  hemodynamic  components, 
including  changes  in  level  of  hemoglobin  oxygenation  (i.e.,  deoxy-  and 


oxy-hemoglobin),  cerebral  blood  flow,  and  the  rate  of  cerebral  oxygen 
metabolism,  for  a  detailed  understanding  of  the  neurovascular  coupling 
phenomenon,  both  during  and  in  between  seizures(Yao  et  al.,  2009; 
Yuan  et  al.,  2007). 

The  present  work  represents  a  major  technological  and  scientific 
improvement  in  epilepsy.  Previous  attempts  to  track  from  moment 
to  moment  populations  of  neurons  participating  in  an  epileptiform 
event  had  not  been  possible  with  any  technique  or  tool  because  of 
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Fig.  6.  Changes  of  blood  vessel  diameter  during  interictal  onset,  (a.)  PA  image  of  the  rat 
cortex  vasculature  with  the  intact  scalp  and  skull.  The  red  arrow  indicates  the 
microvascularity  (MV)  along  the  yellow  dashed  line  direction  for  a  cross-sectional 
scanning  using  a  50  MHz  ultrasound  transducer,  (b.)  Typical  photoacoustic  signals 
from  the  targeted  vessel  at  different  times,  (c.)  EEC  recordings  show  the  interictal 
spikes  and  (d.)  Change  of  the  vessel  size  captured  by  PAT.  Here  was  a  clear  correlation 
between  the  interictal  spikes  and  the  changes  of  the  vessel  size.  Error  bars  (±s.d.) 
were  calculated  from  20  consecutive  spikes.  (Eor  interpretation  of  the  references  to 
color  in  this  figure  legend,  the  reader  is  referred  to  the  web  of  this  article.) 


Spatial  and/or  temporal  sampling  limitations.  The  high  spatial  and 
temporal  sampling  of  the  novel  PAT  system  allowed  for  the  first 
time  the  complete  mapping  of  an  epileptiform  event  in  vivo.  Another 
major  improvement  is  that  this  is  the  first  report  of  mapping  an  epi¬ 
leptiform  event  at  depths  well  below  the  cortex.  In  terms  of  impact, 
this  is  the  first  demonstration  of  emerging  optical  mapping  tool  in 
the  surgical  evaluation  of  focal  cortical  epilepsy,  since  accurate  local¬ 
ization  of  the  epileptic  focus,  propagation  paths,  and  subcortical  net¬ 
works  critically  depends  on  the  precise  localization  of  epileptogenic 
neurons  and  networks.  The  use  of  our  PAT  system,  experimental 


paradigm,  results,  and  analyses  advances  our  understanding  of  epi¬ 
lepsy  and  seizure  temporal  spatial  properties  relative  to  previous  at¬ 
tempts.  The  noninvasive  yet  whole  surface  and  depth  capabilities  of 
the  PAT  system  allowed  for  the  first  time  to  actually  see  what  is  hap¬ 
pening  during  ictogenesis  in  terms  of  seizure  onset  and  spread.  The 
challenge  of  mapping  focal  epilepsy  stems  from  the  observation  that 
the  pathology  associated  with  focal  epilepsy  is  often  distributed 
across  a  number  of  brain  sites.  Seizures  in  animal  models  and  in 
people  often  have  a  multifocal  or  broadly  synchronized  onset.  The  im¬ 
plication  of  these  observations  is  that  each  of  the  sites  could  act  inde¬ 
pendently  to  initiate  a  seizure  or,  potentially,  to  drive  another  site 
into  a  seizure.  The  current  study  lends  support  to  the  theory  that  sei¬ 
zure  onset  and  spread  involves  a  rich  interplay  between  multiple  cor¬ 
tical  and  subcortical  foci  during  the  onset  and  spread  of  focal  epilepsy. 
The  findings  are  timely  in  the  sense  that  the  neuroscience  community 
is  questioning  the  long-held  dogma  that  seizure  onset  involves  some 
sort  of  single  focus  or  epicenter  in  favor  of  the  emerging  thought  that 
seizures  instead  involve  multiple  cortical  and  subcortical  foci. 
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Abstract 

Study  of  neuro-hemodynamic  activities  in  freely  moving  small  animals  (rather  than  anesthetized) 
provides  more  realistic  and  accurate  information  about  diseases,  and  hence  better  understanding. 
While  it  has  been  proved  that  hemodynamic  changes  are  closely  related  to  epileptic  seizures, 
methods  for  detection  in  freely  moving  (small)  animals  are  limited.  In  this  work,  we  integrated 
photoacoustic  sensor  and  EEG  system  into  a  small  device  that  can  be  attached  on  the  head  of 
awaken  rats,  and  for  the  first  time  it  shows  the  realization  of  long  time  and  simultaneous 
monitoring  of  photoacoustic  and  EEG  signals  in  pentylenetetrazol  (PTZ)  induced  seizure  on  freely 
moving  small  animals  i.e.,  rats.  Results  showed  that  both  the  neural  and  vascular  responses  to 
seizure  in  freely  moving  rats  have  characteristics,  which  are  observed  to  be  different  and  more 
diverse  from  that  of  anesthetized  rats,  and  this  calls  for  more  detailed  study  in  future.  This 
technology  also  promises  for  other  hemodynamic  related  research  study  in  freely  moving  small 
animals. 

Introduction 

Epilepsy  is  the  most  common  brain  disease  that  affects  1%  of  the  world  population  [1].  Epileptic 
seizures  are  caused  by  synchronous,  rhythmic  firing  of  a  population  of  neurons  in  the  brain  that 
involves  enormous  increase  in  metabolic  rate  of  oxygen  and  thus  hemodynamic  parameters  such 
as  cerebral  blood  flow  (CBF),  cerebral  blood  volume  (CBV)  and  oxy-saturation  are  greatly 
changed  [2].  Abetter  understanding  of  coupling  between  the  neural  and  vascular  systems  is  crucial 
to  diagnosis,  prediction  and  treatment  for  this  disease. 

Currently  available  neural  imaging  modalities  such  as  functional  magnetic  resonance  imaging 
(fMRI),  positron  emission  tomography  (PET),  single-photon  emission  computed  tomography 
(SPECT),  have  already  been  employed  in  the  study  of  neurovascular  coupling  in  epilepsy  both  in 
human  and  animal  models  [2,  3].  However,  one  remaining  problem  for  these  technologies  is  that 
they  are  expensive  and  bulk  in  size,  so  the  objects  must  be  anesthetized  or  keep  still  during  the 
entire  recording  process.  For  recently  developed  optical  mapping  methods  such  as  optical 
recording  of  intrinsic  signals  (ORIS),  optical  coherence  tomography  (OCT)  and  voltage  sensitive 
dye,  pre-operation  surgeries  like  skull  opening  are  generally  required  due  to  their  limited  imaging 
depth  [4,  5]  and  thus  they  permit  only  anesthetized  animals. 

Problem  associated  with  anesthetized  animals  is:  the  anesthetization  suppresses  the  activities  of 
neurons  and  cardiovascular  system,  so  the  neural  and  vascular  responses  to  seizure  are  altered  [6]. 


As  a  consequence,  awakening  and  freely  moving  animals  with  vascular  response  uninterrupted  by 
extrinsic  factors  is  ideal  and  preferred  model  for  epilepsy  study.  Conventional  EEG  or  telemetric 
EEG  system  with  video  monitoring  is  the  gold  standard  for  seizure  detection  and  quantification, 
and  has  already  been  widely  used  on  freely  moving  animals  [7,  8].  Yet  the  compact  devices,  which 
can  be  firmly  fixed  on  the  animals  and  hence  can  reliably  detect  hemodynamic  changes,  are  still 
missing  and  no  studies  have  been  reported. 

Photoacoustic  technologies  have  been  widely  employed  in  vascular  network  visualization  of  small 
animal  brain,  tumor  detection,  and  in  early  osteoarthritis  diagnosis  [9-12].  In  photoacoustic  effect 
when  a  short  pulse  of  laser  light  is  absorbed  by  tissue,  wideband  ultrasound  signal  that  conveys 
the  position  and  absorbed  energy  of  the  tissue  is  generated  due  to  transient  thermoelastic 
expansion.  In  this  way,  it  compromises  the  merits  of  high  spatial  resolution  and  imaging  depth 
from  ultrasound  imaging,  and  high  intrinsic  optical  contrast  from  optical  imaging.  It  has  been  used 
for  neurovascular  coupling  study  of  epilepsy  with  anesthetized  animals  [9,  10].  However,  most  of 
the  photoacoustic  techniques  focus  on  imaging  that  permits  only  the  anesthetized  animals.  Besides, 
the  imaging  speed  of  photoacoustic  systems  is  often  limited  by  the  laser  repetition  rate.  On  other 
hand,  its  potential  for  real-time  monitoring  as  a  sensor  in  simple  and  compact  design  that  can  be 
used  on  freely  moving  animals  is  not  fully  explored. 

The  purpose  of  this  work  is  to  develop  an  interface  that  consists  of  an  EEG  system  and  a 
photoacoustic  sensor,  and  to  study  the  interaction  between  the  EEG  signal  and  the  hemodynamic 
signal  of  superior  sagittal  sinus  (SSS)  in  pentylenetetrazol  (PTZ)  induced  seizure  on  freely  moving 
rats.  PTZ  has  been  employed  to  induce  generalized  seizure  in  rats  to  study  seizure  phenomenon, 
where  hemodynamic  changes  have  reportedly  been  found  all  over  the  brain.  As  one  of  the  most 
distinguishable  as  well  as  main  blood  vessels  in  brain,  SSS  allows  blood  to  drain  from  lateral 
aspects  of  anterior  cerebral  hemisphere  to  confluence  of  sinuses  and  return  to  venous  circulation, 
covering  large  region  of  brain.  Thus  the  hemodynamic  signal  of  SSS  is  closely  related  to  the  status 
of  the  PTZ-induced  generalized  seizure. 

Through  this  dual-modality,  monitoring  epilepsy  in  freely  moving  animals,  one  can  get  an 
in-depth  knowledge  of  responses  from  vascular  and  neural  systems,  and  hence  better 
understanding  of  mechanisms  of  epilepsy.  This  technology  may  have  greater  utility  in  applications 
such  as  drug  development  and  therapy  assessment  in  epilepsy,  and  other  hemodynamic  related 
studies. 

Material  and  methods 

(1)  System  overview 

The  entire  system  can  be  divided  into  four  main  parts  as:  photoacoustic  system,  EEG  system, 
video  camera,  and  small  interface  mounted  on  rat  head,  as  shown  in  Fig.  1.  Pulse  light  of 
wavelength,  1064nm  from  a  Q-switch  Nd:YAG  laser  operating  at  lOHz  was  coupled  into  a  0.4mm 
core  multimode  fiber  patch  cable  using  a  convex  lens  and  delivered  to  the  rat  brain.  The  transient 
photoacoustic  signal  generated  by  the  short  laser  pulse  was  detected  by  a  pair  of  transducers  with 
central  frequency  5MHz.  The  signal  was  amplified  using  a  26dB  PCB  preamplifier  and  a  17dB 
commercial  preamplifier  (AH- 11 00,  ONDA  Co.),  and  finally  digitalized  by  a  PCI  DAQ  card  (NI 


PCI-5154)  with  sampling  rate  50MHz.  The  intracranial  EEG  signal  was  collected  through  a  pair  of 
small  screws  drilled  in  the  skull  of  rat.  The  depth  of  the  screws  through  the  skull  was  about  1mm. 
The  EEG  signal  was  amplified  (RA16PA,  TDtucker-Davis  Tech.),  and  collected  (RZ5  Bioamp 
Processor,  TDtucker-Davis  Tech.)  and  stored  into  a  computer  for  processing  later.  The  entire 
system  was  synchronized  through  triggering  with  Nd:YAG  laser.  Continuous  video  recording  was 
achieved  with  a  video  camera  (we  used  iPhone  4S).  The  transducers,  multimode  fiber  and  EEG 
connecting  wires  are  detachable  from  the  rat  head  when  the  experiment  is  over.  All  the  data 
including  the  photoacoustic  signal,  EEG  signal  and  video  are  jointly  analyzed  later  for  study  of 
epilepsy. 


Q-switch  Nd:YAG  laser 


Fig.  1:  Schematic  of  experimental  system.  The  integration  of  photoacoustic  system,  EEG  system,  video 
camera,  and  small  interface  mounted  on  the  rat  head  is  depicted. 


(2)  Photoacoustic  sensor  design 

Fig.  2(a)  depicts  the  schematic  of  photoacoustic  sensor  being  mounted  on  a  rat  head.  Core 
multimode  fiber  patch  cable  (0.4mm)  for  delivery  of  light  was  connected  to  an  optical  cannula 
(CFMC14L02,  Thorlabs  Inc.)  with  a  ceramic  mating  sleeve  (see  Fig.  2(b)).  The  cannula  was 
implanted  right  above  SSS,  with  a  tip  (2mm)  reaching  into  brain  through  a  small  hole  in  skull.  In 
rat  experiment,  the  optical  energy  coming  out  of  the  multimode  fiber  was  about  6mJ/pulse  and  that 
from  the  connula  tip  into  the  brain  was  about  4mJ/pulse.  The  cannula  was  fixed  with  dental 
cement,  and  then  the  skull  surface  was  sealed  with  a  thin  layer  of  super  glue.  The  two  transducers 
were  tilted  and  placed  with  their  axes  pointing  to  SSS.  They  were  supported  by  a  3D  printed 
VeroWhitePlus  interface  and  fastened  with  small  screws  and  thermal  glue.  The  space  between  the 
skull  and  the  transducers  was  filled  with  ultrasound  gel  for  proper  acoustic  coupling.  Fig.  2(c) 
shows  a  photograph  of  the  above-mentioned  design  being  mounted  on  rat  head  (after  the  surgery). 


(3)  Surgery 

Four  Sprague  Dawley  rats,  whose  weights  were  around  400g,  were  employed  in  the  experiments. 
SSS  was  about  2mm  below  the  skull.  Fig.  3  shows  the  top  view  of  a  rat  skull  during  surgery.  The 
positions  of  the  optical  cannula  and  the  two  EEG  electrodes  are  noted  by  arrows.  The  surgery  was 


carried  out  as  follows.  Firstly,  rat  was  anesthetized  with  isoflurane  and  hair  on  the  head  was 
removed,  and  subsequently,  the  scalp  was  opened.  Next,  a  hole  of  diameter  0.5mm  was  carefully 
drilled  in  the  skull,  and  cannula  was  inserted  into  the  hole  and  fixed  with  dental  cement.  The 
surface  of  the  skull  was  cleaned  and  applied  with  antibiotic,  and  then  sealed  with  a  thin  layer  of 
super  glue  to  prevent  from  infection.  After  that,  two  small  screws  were  drilled  at  forehead  and  they 
served  as  two  EEG  electrodes.  Finally,  3D  Printed  VeroWhitePlus  interface  was  fixed  onto  the  rat 
skull  with  dental  cement,  and  the  two  EEG  screws  were  also  embedded  in  dental  cement,  leaving 
two  wires  tied  on  the  screws  stretching  out  for  later  connecting  to  recording  leads.  After  four  days 
(for  complete  recovery),  rat  was  adopted  for  (epilepsy)  experiment.  All  the  procedures  were 
carried  out  in  accordance  with  an  approved  University  of  Florida  lACUC  protocol. 
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Fig.  2:  (a)  Schematic  of  the  photoacoustic  sensor  in  3D.  Laser  is  delivered  to  SSS  through  an  optical  cannula, 
which  is  implanted  on  rat  skull.  Then  a  pair  of  transducers  receives  the  generated  ultrasound  signals.  The  diameter 
of  the  receiving  transducers  is  approximately  5.5mm,  and  they  are  supported  by  a  3D  printed  VeroWhitePlus 
interface,  (b)  Photograph  of  optical  cannula,  mating  sleeve  and  multimode  fiber  patch  cable.  The  optical  cannula 
has  a  diameter  of  2.5mm,  with  a  piece  of  core  multimode  fiber  (0.4mm)  in  its  center,  (c)  Photograph  of  rat  head, 
with  all  detectors  being  connected,  after  surgery.  The  EEG  electrodes  are  implanted  in  forehead  of  the  rat,  as 
indicated. 


Fig.  3:  Photograph  of  a  rat  skull  while  in  surgery,  with  the  positions  of  optieal  eannula  and  EEG  eleetrodes  as 
indieated  by  arrows. 


(4)  Data  collection  and  processing 

Phantom  experiments  were  carried  out  to  validate  the  photoacoustic  signal  from  SSS.  A  tube  with 
inner  diameter  0.2mm  was  filled  with  ink,  and  embedded  in  a  phantom  (we  employed  agar 
phantom)  at  a  depth  of  2mm.  Then  VeroWhitePlus  interface  along  with  optical  cannula  and 
transducers  was  placed  onto  the  phantom  to  record  the  photoacoustic  signal  from  the  tube.  By 
comparing  the  resulting  photoacoustic  signal  with  that  of  rat  experiments,  the  photoacoustic  signal 
of  SSS  can  be  verified  (amplitude  vs.  time  was  calculated).  Because  the  optical  absorption  at 
wavelength  1064nm  is  dominantly  sensitive  to  oxy-hemoglobin  (Hb02)  as  compared  to  that  of 
deoxy-hemoglobin  (HbR),  the  photoacoustic  signal  at  this  wavelength  can  be  taken  as  the  volume 
change  of  Hb02  in  SSS  [13].  For  each  experiment,  the  signals  from  the  two  transducers  were 
compared,  and  the  better  one  was  selected  for  recording. 


There  were  two  strong  sources  of  noise  in  the  collected  EEG  data.  One  is  from  the  alternating 
current  supply  at  60Hz  and  other  is  from  Q-switch  of  the  laser  at  lOHz  (Fig.  4(b),  blue  line).  The 
EEG  data  was  firstly  Fourier  transformed  and  the  alternating  current  noise  of  60Hz  was  eliminated 
by  filtering  out  the  frequency  band  in  the  range  58Hz  to  62Hz,  and  the  Q-switch  noise  was 
removed  by  filtering  out  the  frequency  band  of  10  ±  0.5Hz,  20  ±  0.5Hz,  30  ±  0.5Hz. . .  and  so  forth 
(Fig.  4(a)).  Then  the  filtered  EEG  data  was  inverse  Fourier  transferred  back  to  obtain  the  time 
resolved  data  (red  line  as  shown  in  Fig.  4(b)). 


Fig.  4:  EEG  data  de-noising  in  frequency  domain,  (a)  Frequency  power  chart  of  EEG  signal  (before  de-noise  (red) 
and  after  de-noise  (blue)),  (b)  EEG  signal  before  (blue)  and  after  (red)  de-noise.  The  noises  are  indicated  with 
arrow  and  green  encircles:  A  (the  Q-switch  noise  (lOHz)  from  the  laser)  and  B  (the  alternating  power  supply  noise 
(60Hz)). 


(5)  Rat  experiments 

All  the  rat  experiments  were  carried  out  between  2:00pm  and  5:00pm.  The  rats  were  first  placed  in 
an  acrylic  cage  and  acclimatized  for  at  least  one  hour  before  the  experiment.  Then  the  rats  were 
anaesthetized,  and  the  EEG  leads,  the  multimode  fiber,  and  the  transducers  were  connected  to  the 
head  of  rats  for  recording. 

For  free  moving  rat,  the  seizure  experiment  was  initiated  around  lOmin  after  the  awake  of  rat,  in 


which  seizures  were  generated  by  injecting  PTZ  of  30-75mg/Kg  into  abdomen  [14].  Two  kinds  of 
control  tests  were  performed.  In  first  case,  sterile  saline  of  same  volume  was  injected,  and  in  other 
case,  rats  were  continuously  anesthetized  with  a  PTZ  dose  of  120mg/Kg.  The  concentration  for 
PTZ  was  fixed  to  50mg/mL.  The  duration  of  recording  is  30mins  for  each  experiment,  and  the 
injection  of  drug  was  induced  around  5mins.  Behavioral  seizures  of  all  rats  were  scored  using 
5-graded  Racine  Score  system  [15],  and  the  score  was  based  on  the  severest  behavior  of  rats.  A 
total  number  of  22  experiments  (n=22)  were  carried  out  using  4  rats,  as  shown  in  Table  1 .  The 
time  interval  between  two  experiments  for  the  same  rat  was  at  least  3  days.  All  the  results  were 
confirmed  by  experts  in  neurology. 

Table  1:  Experiment  design.  The  number  in  braekets  indieates  the  number  of  rats  that  were  sueeessfully  video 
reeorded  through  SOmins.  The  star  mark  indieates  that  the  rat  died  after  experiment. 


case  number 

rati 

rat2 

rats 

rat4 

anaesthetized 

1* 

1* 

1 

0 

30mg/Kg  PTZ 

2(2) 

1(1) 

0 

1(1) 

50mg/Kg  PTZ 

2 

4(1) 

3(2) 

0 

75mg/Kg  PTZ 

0 

0 

0 

1 

control 

2 

1 

2 

0 

Results 

(1)  Photoacoustic  signal  verification 

Photograph  of  a  phantom  (employed  in  experiment)  is  shown  in  Fig.  5(a)  and  photoacoustic 
signals  for  phantom  experiment  (blue  line)  and  rat  experiment  (red  line)  are  shown  in  Fig.  5(b). 
Both  the  curves  show  strong  and  sharp  bipolar  signals  at  the  instant  around  3us  with  clear 
background,  which  is  equivalent  to  about  5.3mm  in  tissue.  We  conclude  that  the  photoacoustic 
signal  from  SSS  is  about  3us  in  rat  experiment. 


Fig.  5  Phantom  experiment  for  validating  photoaeoustie  signal  of  SSS.  (a)  Photograph  of  phantom,  (b) 
Photoaeoustie  signals  in  phantom  experiment  (blue),  and  in  rat  experiment  (red).  The  signal  for  SSS  is  at  instant 
about  3.5us. 


(2)  Dose  and  behavior 

In  all  unanesthetized  rat  experiments,  the  rat  was  put  into  an  uncovered  acrylic  cage,  as  shown  in 
Fig.  6.  The  dose  of  PTZ  was  carefully  administrated  ranging  from  30-75mg/Kg.  One  rat  died 


shortly  after  a  dose  of  75mg/Kg,  while  all  the  rats  injected  with  50mg/Kg  or  under  survived. 
During  the  seizure  process,  the  rats  showed  behavior  of  facial  clonus,  head  nodding,  forelimb 
clonus,  rearing  and  falling  (loss  of  postural  control;  See  supplemental  video  S.  1).  The  video  was 
recorded  about  0.5m  in  front  of  the  cage  with  iPhone  4S.  In  some  experiments,  the  rats  in  the  cage 
were  well  trapped.  In  this  case,  the  rats  only  turned  around  the  two  ends  of  the  cage,  and  then  the 
operator  needed  to  turn  the  cage  in  opposite  direction  to  prevent  the  multimode  fiber  from  twisting. 
However,  in  other  experiments,  the  rats  jumped  out  of  the  box,  or  tom  off  the  optical  fiber  and/or 
transducers.  Thus,  only  33%  (3  out  of  10)  of  the  freely  moving  rat  experiments  with  more  than 
50mg/Kg  PTZ  can  be  fully  recorded  all  through  the  experiment  (Table  1).  No  seizure  was 
observed  in  the  rats  injected  with  sterile  saline.  In  general,  the  severity  of  seizure  tended  to  be 
higher  with  a  higher  dose  of  PTZ  injection  (see  Fig.  7).  For  anesthetized  rat  experiments,  the  rats 
were  continuously  anesthetized  with  isoflurane  breathing  circuit,  and  the  body  temperature  of  the 
rats  was  maintained  at  37  V  with  a  thermo  pad. 


Fig.  6:  Awake  and  freely  moving  rat  in  aerylie  eage  for  experiments. 


Racine  score  Dill  ^  IV  [Z]  V 


Fig.  7:  Distribution  of  Raeine  seore  (in  %)  in  freely  moving  rat  seizure  experiments  for  eaeh  PTZ  dose. 

(3)  Simultaneous  EEG  and  hemodynamic  response  study 

Fig.  8  shows  the  experiment  results  of  three  different  rat  models:  anaesthetized  seizure  model  (Fig. 
8(a)),  control  model  (Fig.  8(c)),  and  freely  moving  seizure  model  (Fig.  8(b)  using  rat2  and  Fig.  8(d) 
using  rat3)  respectively.  All  the  injections  were  administrated  at  the  instant  5min.  By  comparing 
these  curves,  we  noted  that  these  three  different  rat  models  showed  distinguishable  features  for 
both  Hb02  (red)  and  EEG  (blue)  signals. 


For  anesthetized  rats  (with  120mg/Kg  PTZ  injection,  Fig.  8(a)),  there  are  few  sharply  isolated 
waves  in  the  EEG  signals  before  the  administration  of  PTZ,  and  at  the  same  time  some 
fluctuations  in  Hb02  are  observed.  The  PTZ  injection  causes  some  high  amplitude  spikes  in  the 
EEG  signal  at  5min.  After  seizure  begins,  the  EEG  signal  slowly  increases  in  amplitude  and 
frequency,  and  in  the  mean  while,  the  Hb02  signal  first  shows  a  sudden  drop,  and  then  gradually 
increases.  Sooner  after,  some  hypoxia  gaps  appear  in  the  Hb02  signal,  which  becomes  wider  and 
deep  as  seizure  evolves.  It’s  noted  that  two  anesthetized  rats  were  killed  due  to  overdose  of  PTZ 
when  they  got  awake  after  taking  isofiurane  breathing  circuit  away  with  completion  experiments. 

For  freely  moving  rats  (with  50mg/Kg  PTZ  injection)  in  Fig.  8(b),  the  EEG  signal  is  active  before 
the  PTZ  injection,  while  Hb02  level  is  relatively  stable.  When  the  rats  were  caught  and  taken  out 
for  injection  at  5min,  the  EEG  signal  activity  becomes  low,  and  the  Hb02  signal  increases 
remarkably  for  a  short  interval  of  time.  The  seizure  onset  is  about  3min  after  the  injection.  At  the 
same  time,  the  Hb02  level  increases  significantly  and  then  gradually  gets  lower,  and  the  EEG 
signal  displays  continuous  high  amplitude  that  shows  bursting  activities.  The  result  in  Fig.  8(d) 
(with  0.5mg/Kg  PTZ  injection)  is  much  different  from  that  in  Fig.  8(b).  The  fluctuation  of  Hb02 
signal  in  Fig.  8(d)  is  much  larger  (all  through  the  experiments)  than  that  in  Fig.  8(b).  The  overall 
activity  of  EEG  signal  is  prominently  increased  due  to  seizure  onset  in  Fig.  8(d),  but  not  in  Fig. 
8(b).  Apparently,  the  results  for  freely  moving  rat  experiments  seemed  more  individual  and  case 
dependent. 

For  the  test  experiment  under  control  (Fig.  8(c)),  in  which  the  rats  were  injected  with  0.4mL  of 
0.9%  sterile  saline,  no  dramatic  change  is  seen  both  in  Hb02  and  EEG  signals  after  the  PTZ 
injection  as  compared  to  the  signals  before  injection.  It’s  also  noted  that  the  injection  itself  doesn’t 
influence  much  either  the  Hb02  signal  or  EEG  signal  as  shown  in  Fig.  8(c),  which  is  much 
different  from  the  results  for  freely  moving  rat  experiments  in  Fig.  8(b)  and  Fig.  8(d).  This  may  be 
due  to  individual  difference  of  rats. 


Fig.  8:  EEG  (blue)  and  Hb02  (red)  signals  for  rat  experiments  of  different  models,  (a)  Anaesthetized  rat  with 
120mg/Kg  PTZ  injeetion.  (e)  Freely  moving  rat  with  0.4mL  injeetion  of  0.9%  sterile  saline,  (b)  and  (d)  Freely 


moving  rats  with  50mg/Kg  PTZ  injection. 


Discussion 

(1)  The  photoacoustic  sensor  design 

Photoacoustic  technologies  have  gained  much  attention  for  their  diverse  applications  in 
biomedical  imaging,  but  it’s  seldom  used  as  a  single  sensor.  In  this  work,  we  aim  to  integrate  a 
photoacoustic  sensor  and  EEG  system  in  a  small  device  for  the  dual  modality  study  of  epilepsy  in 
free  moving  rats.  By  utilizing  ultrasound  as  means  of  detection,  the  penetration  depth  of  PAT  can 
be  much  deeper  than  that  of  pure  optical  imaging  while  keeping  high  spatial  resolution.  Thus 
non-invasive  detection  of  blood  vessels  in  small  rat  brain  can  be  achieved.  However,  in  this  work, 
small  hole  is  required  to  be  drilled  in  rat  skull  because  laser  light  from  multimode  fiber  can’t 
penetrate  scalp  and  skull  of  mature  rat.  Small  rat  can  not  be  employed  here  because  fiber  and 
transducers  can  not  be  firmly  attached  to  the  loose  scalp  of  small  rats.  There  are  other  limitations 
in  this  work  as  well.  For  example,  we  only  used  one  wavelength,  but  it  can  be  easily  adapted  to 
dual  wavelength  measurement  with  the  same  Q-switch  Nd:YAG  laser;  and  multiple  illumination 
and  detection  sites  in  the  brain  can  also  be  realized  too.  All  these  can  be  improved  in  future  work. 

(2)  Dual  modality  study  with  video  monitoring 

When  rat  was  anesthetized,  the  fiber,  transducers,  and  EEG  leads  were  swayed  but  no  notable 
fluctuations  was  observed  in  both  the  photoacoustic  and  the  EEG  signals,  proving  the  reliability  of 
the  system.  Video  monitoring  was  used  to  confirm  the  presence  of  a  motor  component  associated 
with  an  EEG  seizure  and  to  score  the  severity  of  motor  seizures.  We  noticed  that  there  are  some 
artifacts  in  both  the  photoacoustic  signal  and  EEG  signal  before  the  injection  of  drug,  which  are 
absolutely  not  related  to  seizure.  These  artifacts  could  not  be  discriminated  without  video 
monitoring. 

For  simultaneous  monitoring  of  EEG  and  photoacoustic  signals,  the  results  for  anesthetized  rats 
are  generally  similar  in  which  hypoxia  gaps  can  be  easily  observed  in  the  Hb02  signal  in  all  the 
cases.  In  contrast,  results  in  freely  moving  rats  are  more  case  dependent  and  complicated  than  that 
of  anaesthetized  rats.  In  some  experiments  with  freely  moving  rats,  we  also  observed  the  proceeds 
of  hemodynamic  signal  to  the  EEG  signal;  and  EEG  signals  of  different  frequencies  are  decoupled 
sometimes.  However,  in  this  study,  no  statistic-based  analysis  is  derived  due  to  limited  number  of 
cases.  The  EEG  and  hemodynamic  signal  recordings  after  seizure  (which  we  missed  to  carry  out 
in  this  work)  are  also  important  and  interesting  issues  for  epilepsy  study.  As  a  result,  more  detailed 
study  is  required  (in  future)  with  more  cases  and  longer  observation  time.  But  it  is  beyond  the 
scope  of  this  communication. 

Conclusion 

This  study  demonstrates  the  practicality  of  photoacoustic  detection  combined  with  EEG  and  video 
monitoring  of  epilepsy  in  freely  moving  rats.  Rats  with  different  seizure  model  were  tested,  and 
we  found  that  the  EEG  and  hemodynamic  signals  during  seizure  in  freely  moving  rats  are  more 
complicated  than  that  of  anesthetized  rats,  which  requires  more  detailed  research  in  future.  With 
the  detachable  design  of  fiber  and  transducers  in  this  type  of  recording  system,  one  rat  can  be  used 
many  times,  allowing  long-term  studies  of  epilepsy.  In  our  experiments,  rats  survived  more  than 


one  month  until  they  were  killed.  With  appropriate  modifications,  this  technology  can  easily  find 
other  applications  related  to  hemodynamic  study  than  epilepsy,  in  which  freely  moving  rat  model 
is  preferred. 
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Abstract — Neurovascular  coupling  in  epilepsy  is  poorly 
understood,  the  study  of  which  requires  simultaneous 
monitoring  of  hemodynamic  changes  and  neural  activity  in  the 
brain.  Here  we  for  the  first  time  present  a  combined  real-time 
3D  Photoacoustic  tomography  (PAT)  and 
electro-physiology/electroencephalography  (EEG)  system  for 
the  study  of  neurovascular  coupling  in  epilepsy,  whose  ability 
was  demonstrated  with  a  pentylenetetrazol  (PTZ)  induced 
generalized  seizure  model  in  rats.  Two  groups  of  experiments 
were  carried  out  with  different  wavelengths  to  detect  the 
changes  of  oxy-hemoglobin  (Hb02)  and  deoxy-hemoglobin 
(HbR)  signals  in  the  rat  brain.  We  extracted  the  average  PAT 
signals  of  the  superior  sagittal  sinus  (SSS),  and  compared  them 
with  the  EEG  signal.  Results  showed  that  the  seizure  process 
can  be  divided  into  three  stages.  A  ‘dip’  lasting  for  1-2  minutes 
in  the  first  stage  and  the  following  hyperfusion  in  the  second 
stage  were  observed.  The  oxy-hemoglobin  (Hb02)  signal  and 
the  deoxy-hemoglobin  (HbR)  signal  were  generally  negatively 
correlated.  The  change  of  blood  flow  was  also  estimated.  All  the 
acquired  results  here  were  in  accordance  with  other  published 
results.  Compared  to  other  existing  functional  neuroimaging 
tools,  the  proposed  method  here  enables  reliable  tracking  of 
hemodynamic  signal  with  both  high  spatial  and  temporal 
resolution  in  3D,  so  it’s  more  suitable  for  neurovascular 
coupling  study  of  epilepsy. 

1.  Introduction 

Neurovascular  coupling  refers  to  relationship  between 
neuronal  activity  and  hemodynamic  related  changes  in  the 
brain,  which  are  closely  coupled  in  time,  spatial,  and 
amplitude.  In  epileptic  seizures,  the  synchronous  excessive 
neuron  discharges  caused  an  enormous  increase  in  the 
metabolic  rate  of  oxygen,  which  will  place  supra-normal 
demands  on  the  brain’s  auto-regulatory  mechanisms,  thus  the 
normal  neurovascular  coupling  can  be  greatly  changed  [1,2]. 
In  order  to  precisely  characterize  and  understand  the 
mechanisms  of  neurovascular  coupling  during  seizures, 
which  is  still  poorly  understood,  and  potentially  benefit  the 
predication  and  diagnosis  of  epilepsy,  simultaneous 
recording  of  hemodynamic  related  signals  and  neuronal 
activities  is  required. 

Although  some  investigators  have  tried  new  technologies 
such  as  membrane  voltage  sensitive  dye  to  image  the  animal 


cortex  [3],  EEG  still  serves  as  the  most  common  and  classic 
way  for  neural  activity  recording  [4].  The  biggest  problem 
with  current  practice  of  the  neurovascular  coupling  study  is 
the  lack  of  an  effective  way  to  read  out  the  hemodynamic 
changes  in  the  whole  brain  with  high  spatiotemporal 
resolution.  The  disadvantages  with  currently  available 
neuroimaging  technologies  that  detect  the  hemodynamic 
related  signals  are  obvious.  Traditional  clinical  modalities 
such  as  functional  magnetic  resonance  imaging  (fMRI), 
positron  emission  tomography  (PET),  and  single-photon 
emission  computed  tomography  (SPECT)  suffer  from  their 
limited  spatial  and  temporal  resolutions  [5,  6];  while  more 
recently  optical  mapping  methods  such  as  recording  of 
intrinsic  signals  (ORIS)  can  only  see  superficial  cortical 
changes  due  to  the  limited  penetration  depth  [7].  Thus  there  is 
a  growing  interest  in  developing  new  techniques  for 
neurovascular  coupling  study,  a  good  understanding  of  which 
will  in  return  help  the  interpretation  of  the  acquired  results  of 
these  modalities. 

Photoacoustic  tomography  (PAT)  is  a  hybrid  method  that 
is  capable  of  imaging  optical  absorption  of  tissue  through  the 
detection  of  ultrasound  waves  generated  by  a  short  laser  pulse 
due  to  transient  thermoelastic  expansion  [8-1 1].  It  enables  the 
label-free  vascular  imaging  with  a  high  signal-to-noise  ratio 
(SNR)  by  utilizing  the  optical  absorption  of  hemoglobin  as 
contrast,  the  primary  carrier  of  oxygen  in  the  blood.  Besides, 
compared  with  pure  optical  imaging,  it  uses  less  scattered 
ultrasound  as  its  probing  signal  and  thus  offers  a  much  deeper 
penetration  and  offers  3D  imaging  ability  with  high  spatial 
resolution.  With  high  repetition  laser  and  multi-wavelength 
method,  hemodynamic  parameters  such  as  oxy-saturation, 
total  blood  volume,  and  even  blood  flow  can  also  be  detected 
in  real-time.  With  all  these  merits,  PAT  is  more  suitable  for 
neurovascular  study. 

PAT  can  be  implemented  in  many  ways,  but  most  of  these 
systems  are  based  on  2D.  To  get  3D  images,  mechanical 
scanning  along  one  direction  is  needed  [8,  9],  and  the  frame 
rate  is  limited  by  the  repetition  rate  of  the  laser  used.  Besides, 
some  of  these  systems  will  have  a  limited  imaging  depth  or 
scanning  area  [9],  while  most  seizure  involves  rapid  changes 
of  widespread  cortical  networks  in  3D  [12].  The  most 
effective  way  for  real-time  3D  imaging  is  to  use  2D 
transducer  array  with  multi-channel  acquisition  system. 
Although  2D  planar  array  is  easy  for  fabrication,  it  has  a  low 
transverse  resolution  due  to  the  limited  aperture,  and  the 
features  with  high  aspect  ratio  or  with  orientations  oblique  to 
the  transducer  surface  will  suffer  from  distortion  [13]. 
Comparatively,  a  spherical  array  can  offer  more  complete 
angular  views  of  the  object,  providing  both  high  resolution 
and  accurate  feature  definition  regardless  of  shape  or  location 
of  the  object. 

Small  animal  brain  imaging  is  one  of  the  major  directions 
of  the  PAT,  but  most  works  there  are  mainly  focused  on 
vascular  structural  imaging,  tracking  single  blood  vessel 


signals,  or  brain  functional  mapping  [9],  and  research  in 
epilepsy  is  relatively  rare.  PAT  in  2D  with  concurrent  EEG 
for  epilepsy  study  has  been  reported,  but  the  conclusions  in 
conflict  with  the  existing  knowledge  may  be  due  to  wrong 
interpretation  of  the  data  [14].  Previously,  we  for  the  first 
time  reported  a  sparse  spherical  array  based  PAT  system 
operating  at  a  frame  rate  of  3.3f/s  for  focused  seizure 
monitoring  but  without  concurrent  EEG  [10].  Here  we 
modified  the  system  by  incorporating  EEG  system  into  it  for 
the  study  of  neurovascular  coupling  of  pentylenetetrazol 
(PTZ)  induced  generalized  seizure  in  rats.  Two  groups  of 
experiments  were  conducted  at  two  different  wavelengths  to 
test  the  combined  system,  and  the  acquired  PAT  and  EEG 
datas  were  jointly  analyzed. 


animal  holder  using  glue  tapes,  so  as  to  avoid  noise  induced 
by  any  possible  unnecessary  movement  of  the  rat  or 
experiment  operators.  The  PAT  system  gave  an  output 
synchronized  signal  to  the  EEG  system,  so  that  the  two 
systems  were  well  synchronized. 


II.  MATERIALS  AND  METHODS 

A.  System  setup 

The  real-time  3D  photoacoustic  system  has  already  been 
reported,  and  the  main  feature  of  this  system  is  the  sparse 
spherical  array  which  is  consisted  of  192  discrete  transducers, 
as  shown  in  Fig.  1.  The  transducers  were  mounted  on  a 
custom  fabricated  white  ABS  spherical  interface  with  a 
distance  of  71.6mm  from  the  center  of  the  array,  and  formed 
into  seven  layers.  Each  transducer  had  a  central  frequency  of 
5MHz  and  a  reception  bandwidth  of  greater  than  80%, 
resulting  in  an  almost  isotropic  spatial  resolution  about 
0.2mm.  The  active  area  of  the  transducer  was  3mm  in 
diameter,  and  the  angular  acceptance  was  about  15  degree. 
The  signal  from  the  192  transducers  was  amplified  by  16 
homemade  preamplifier  boards  and  coupled  into  a  64  channel 
parallel  data  acquisition  system  with  3:1  multiplexing,  as 
indicated  with  three  different  colors. 

Rats  were  elevated  to  the  center  of  the  spherical  interface 
through  a  chamber  fixed  at  the  tank  bottom,  whose  top  was 
about  1 5mm  beneath  the  interface  center,  and  a  transparent 
plastic  wrap  was  used  to  cover  the  chamber  top.  Ultrasound 
gel  was  used  to  couple  the  sound  from  the  animal  to  the 
plastic  wrap.  The  ultrasonic  transducer  array  with  the 
spherical  interface,  placed  in  the  water  bath,  was  therefore 
acoustically  coupled  to  the  mouse  tissues.  The  laser  light  was 
delivered  through  a  concave  lens  giving  a  homogenous 
illumination  on  the  head  of  the  rat.  In  this  study,  two 
wavelengths  were  used,  1064nm  and  710nm.  The  1064nm 
light  was  from  a  pulsed  Nd:YAG  (5~20ns  pulse  width),  and 
the  710nm  light  was  from  a  Q-switched  Nd:YAG  pumped 
Ti:  sapphire  laser.  The  final  laser  intensity  on  the  rat  head  was 
about  50mJ/cm^  for  1064nm,  and  4mJ/cm^  for  710nm,  all 
below  the  ANSI  exposure  limits.  Both  the  two  lasers  operate 
at  lOHz,  so  that  the  system  could  record  one  complete  set  of 
3D  PAT  data  in  0.33s. 

Two  injection  needles,  inserted  into  the  front  cortex  of  rat 
brain  about  3mm  below  the  scalp,  were  served  as  EEG 
electrodes  as  shown  in  Fig.  1(b).  One  electrode  was  picked 
for  reference  and  ground  signal,  while  the  other  electrode  was 
to  record  EEG  signal.  The  EEG  signal  was  amplified 
(RA16PA,  TDtucker-Davis  Tech.)  and  recorded  (RZ5 
Bioamp  Processor,  TDtucker-Davis  Tech.)  with  a  sampling 
rate  about  50kHz.  Both  the  two  electrodes  were  firmly  glued 
on  the  rat  head  and  the  wires  were  tightly  attached  to  the 


Figure  1 .  Experimental  setup  of  the  combined  system,  (a)  Schematic  of  the 
photoacoustic  system;  (b)  Rat  with  two  electrodes  in  the  brain;  and  the  laser 
illumination  is  also  indicated 

B.  Animal  preparation 

Two  groups  of  rats  (n=16  rats,  each  about  35-45g,  with  8 
rats  for  1064nm  and  the  other  8  rats  for  710nm)  were  used. 
The  rats  were  anaesthetized  with  about  0.25ml  Urethane 
(0.25g/mL)  and  got  hair  on  the  head  shaven.  Before  mounting 
the  rats  on  the  homemade  plastic  holder,  a  syringe  containing 
0.04g/ml  PTZ  with  a  scalp  needle  was  inserted  in  the 
abdomen  and  fixed  with  glue.  PTZ  was  used  to  induce 
generalized  seizure  in  rat  brain,  which  has  already  been 
established  [15].  Then  the  two  electrodes  were  inserted,  and 
the  rats  were  elevated  into  the  transducer  array  center  of  the 
photoacoustic  system  for  scanning.  Each  experiment  lasted 
25min,  and  PTZ  was  administrated  with  a  dose  of  0.5mL  at 
5min.  All  rats  were  kept  alive  through  the  experiment.  All  the 
animal  procedures  performed  were  in  accordance  with  the 
approved  University  of  Florida  lACUC  protocols. 

C.  Data  processing 

With  a  frame  rate  of  3.3f/s,  5000  complete  sets  of  PAT 
data  were  recorded.  Before  doing  3D  reconstruction  and 
processing,  2D  PAT  images  were  reconstructed  with  data 
collected  by  the  64  transducers  in  the  lowest  layer  of  the 
array,  which  were  indicated  with  blue  in  Fig.  1.  The  2D 
images  were  carefully  checked  to  avoid  the  possible  motion 
artifact  induced  by  the  movement  of  the  rat  head  during  the 
experiment.  3D  images  were  calculated  using  delay&sum 
algorithm,  which  was  accelerated  by  GPU  parallel  technique. 
Then  the  negative  pixel  values  in  were  set  to  zero.  Each  set  of 
the  reconstructed  3D  images  had  a  voxel  number  of 
201x201x101  representing  a  10xl0x5mm^  volume,  which 
could  be  calculated  in  0.97s.  Amira  (from  Visage  Imaging, 
Inc.)  was  used  to  render  the  3D  images. 

There  were  two  sources  of  noise  in  the  EEG  signal,  a 
strong  sharp  lOHz  noise  from  the  Q-switch  of  the  laser,  and 
the  60Hz  alternating  current  noise.  For  de-noising,  the 
deteriorated  data  by  the  Q-switch  was  replaced  with  adjacent 
data,  and  the  60Hz  alternating  current  noise  was  removed  by 
filtering  out  58-62Hz  in  frequency  domain  (Fig.  2).  Then  the 
EEG  data  was  applied  with  a  low  pass  filter  of  lOOHz,  and 


transformed  into  a  300000  point  signal  for  the  25min  EEG 
recording,  corresponding  to  sampling  rate  of  200Hz. 

For  neurovascular  coupling  study,  we  first  extracted  the 
averaged  PAT  signals  of  the  superior  sagittal  sinus  (SSS)  of 
each  rat  from  the  reconstructed  3D  images  and  compared  it 
with  the  corresponding  EEG  signal.  We  choose  the  PAT 
signals  of  the  SSS  because  it  has  the  highest  SNR  and  is 
highly  distinguishable  in  the  reconstructed  3D  image. 
Because  1064nm  and  710nm  are  mainly  absorbed  by 
oxy-hemoglobin  (Hb02)  and  deoxy-hemoglobin  (HbR),  the 
acquired  results  for  these  two  wavelengths  can  be  taken  as  the 
changes  of  Hb02  and  HbR  respectively  [16].  All  the  acquired 
results  were  examined  by  experts  in  the  neurology. 

III.  RESULTS 


A.  Data  preprocessing 

Representative  time  frames  of  2D  images  through  the 
experiment  are  displayed  in  Fig.  2(a)  with  a  time  interval  of 
Imin,  compared  with  a  rat  brain  photo  taken  after  the 
experiment  (Fig.  2(b)).  Although  the  image  quality  is  not 
excellent  due  to  the  limited  transducer  positions  (only  64 
transducers  are  averaged  once),  the  SSS  and  two  transverse 
sinuses  (TS)  are  clearly  reconstructed,  and  some  other  small 
blood  vessels  are  also  revealed.  In  Fig.  2(a)  it  is  observed  that 
there  is  no  palpable  shift  of  the  rat  brain  which  is  a 
prerequisite  for  the  following  data  processing.  In  our 
experiments,  3  out  of  8  rats  with  I064nm  moved,  but  all  the 
rats  in  the  second  group  remained  still.  We  also  notice  that 
it’s  hard  to  tell  when  the  seizure  begins  just  by  eyes  from 
these  2D  images,  so  quantitative  analyze  is  needed. 


EEG  signals  were  recorded  to  confirm  the  onset,  duration, 
and  morphology  of  the  seizure,  which  is  shown  in  Fig.  3(a). 
The  EEG  signal  is  stable  before  the  drug  injection  at  5min, 
when  isolated  high  amplitude  spikes  can  be  easily  seen  in  the 
EEG  signal.  After  this  spontaneous  fast  discharges  begin  to 
show  and  increase  in  amplitude  in  the  following  2  to  3 
minutes.  The  seizure  onset  comes  about  142  +  45s  after  the 
drug  injection,  which  is  defined  to  be  the  first  signal  that 
exceeds  6  times  the  standard  deviation  (SD)  of  the  baseline 
(the  average  signal  in  the  first  4min).  As  the  seizure  evolves, 
the  firing  frequency  generally  gets  lower,  but  the  EEG  signal 
begins  to  decrease  in  and  keeps  firing  through  the 
experiment.  Fig.  3(b)  shows  the  EEG  signal  in  a  short  time 
window  (blue),  compared  with  the  original  data  before 
preprocessing  (red).  We  can  see  that  both  the  noises  from  the 
Q-switch  and  the  alternating  power  supply  are  removed. 


Figure  2.  Representative  time  frames  of  2D  PAT  images  of  the  rat  brain, 
compared  with  an  anatomical  image,  (a)  Representative  time  frames  of  2D 
PAT  images  of  the  rat  brain  with  1064nm.  The  time  interval  is  Imin,  and 
each  image  covers  10  X  lOmm^  with  201  X  201  pixels,  (b)  Photo  of  rat  brain 
with  hair  and  scalp  being  removed.  The  SSS,  TS,  and  the  reconstruction  area 
is  indicated.  Scale  bar  is  5mm  through  the  figure. 
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Figure  3.  EEG  signal  for  PTZ-induced  seizure,  (a)  Typical  EEG  signal 
after  processing,  (b)  The  original  (blue)  and  processed  (red)  EEG  signal  in  a 
short  time  wiondow.. 


B.  3D  image  reconstruction 

Fig.  4  shows  the  reconstructed  3D  images  from  different 
views.  The  reconstructed  volume  is  10  x  10  x  5mm^  with  a 
0.05mm  voxel  size.  Similar  to  the  2D  images,  the  SSS  and  TS 
can  be  clearly  revealed  both  in  the  3D  images  of  Hb02  and 
HbR,  although  the  SNR  of  HbR  is  not  as  good  as  that  of  Hb02 
due  to  relatively  low  laser  intensity  of  710nm,  as  well  as  the 
low  concentration  of  HbR  in  the  blood  vessels  compared  to 
that  of  Hb02.  It  is  also  noticed  that  in  the  3D  images  of  Hb02, 
we  can  see  the  Ostia  end  of  the  straight  sinus  (SS)  below  the 
SSS  which  can  not  be  seen  in  the  2D  images.  All  these 
demonstrate  the  excellent  positioning  ability  of  our  3D  PAT 
system. 


Figured.  Reconstructed  3D  PAT  images  of  the  rat  brain  with  1064nm 
from  different  views.  SSS,  TS,  and  SS  can  be  seen  in  the  image.  The  image 
domain  is  10  X  10  X  5nim^  with  201  X  201  X  101  pixels.  Colarbar  is 
indicated. 


C.  Neurovascular  coupling  study 

For  each  set  of  PAT  data,  the  average  signal  in  a  ROI  that 
could  enclose  the  SSS  was  extracted  and  taken  as  the  change 
of  the  SSS  (Fig.  5(a)  and  Fig.  5(b)),  which  was  further 
converted  into  percentage  changes  to  the  baseline.  For  all  the 
rat  datas  analyzed  (n=14),  the  EEG  signals  looked  the  same  in 
shape  and  timing  (Fig.  5(c)  and  Fig.  5(d)).  This  means  our 
PTZ  induced  generalized  seizure  modal  is  quite  stable  and  the 
drug  administration  is  well  controlled,  so  that  all  the  rat 
results  are  comparable.  The  two  PAT  curves  in  Fig.  5(c)  and 
Fig.  5(d)  are  also  representatives.  The  PAT  signals  for  HbR  is 
smoothed  using  a  moving  window  of  3  seconds  to  improve 
the  SNR.  The  neurovascular  coupling  of  our  rat  model  can  be 
obtained  by  comparing  these  two  groups  of  PAT  curves  and 
the  corresponding  EEG  signals. 


Figure  5.  Typical  results  for  Hb02  (a  and  c)  and  HbR  (b  and  d) 
respectively,  with  the  ROI  for  the  SSS  indicated  in  the  3D  PAT  images.  The 
reconstruction  domain  for  the  3D  images  (a  and  c)  are  lOxlQxSmm^  with 
201x201x101  pixels,  c  and  d  are  the  EEG  and  PAT  signals  of  SSS  for  a  and 
b  respectively.  The  EEG  signals  (blue  line)  are  in  arbitrary  units,  and  the 
PAT  curves  (red  line)  are  represented  as  the  change  with  respect  to  the 
baseline.  The  whole  seizure  process  after  the  injection  of  PTZ  can  be 
divided  into  three  stages.  The  PAT  signal  for  HbR  has  been  smoothed  using 
a  moving  window  of  3  seconds.  The  time  units  in  c  and  d  is  minutes. 

Both  the  EEG  and  PAT  signals  in  these  two  datas  remain 
quite  stable  before  the  drug  injection.  After  that  the  time 
course  of  the  experiment  can  be  generally  divided  into  three 
stages,  as  indicated.  In  the  first  stage,  the  EEG  signal 
increases  in  amplitude  as  seizure  initiates.  At  the  same  time, 
there’s  a  prominent  sudden  decrease  in  the  Hb02  signal  while 
the  HbR  signal  just  slightly  increases.  This  indicates  that  for 
this  period  of  time,  the  increase  of  blood  supply  was 
inadequate  for  the  metabolism  of  neurons.  In  the  second  step, 
we  can  see  that  the  EEG  signals  continue  to  increase  as  the 
seizure  evolves.  However,  the  Hb02  signal  begins  to  increase 
and  the  HbR  begins  to  decrease  at  the  same  time.  Thus  we 
conclude  that  the  blood  supply  is  greatly  increased  due  to 
dilation  of  arteries  and  increasing  of  blood  flow.  In  the  third 
step,  although  the  EEG  signal  gradually  decreases  in  firing 
frequency,  the  Hb02  signal  is  still  decreased,  while  the  HbR 
signal  increases.  This  means  the  blood  flow  is  decreased. 

Through  the  whole  seizure  process,  we  see  the 
coincidence  that  the  two  PAT  signals  are  well  negatively 
correlated.  When  the  Hb02  reaches  its  peak  the  HbR  just 
drops  to  its  minimum,  or  otherwise.  To  make  it  more 
convincing,  we  calculate  the  average  changes  of  Hb02  and 


HbR  across  all  the  rat  experiments,  as  shown  in  Fig.  6.  The 
two  curves  there  are  achieved  by  first  normalizing  the  data  in 
a  single  experiment  to  the  baseline,  and  then  averaging  the 
data  of  all  rats  in  the  same  group.  The  final  data  are  further 
averaged  every  30secs,  with  the  corresponding  SD  indicated. 
The  correlation  coefficient  of  these  two  signals  is  calculated 
to  be  -0.705. 


Figure  6.  Full  timecourses  of  mean  percent  changes  of  the  Hb02  (a)  and 
HbR  (b)  signals,  which  are  negatively  correlated  with  each  other.  The  SD 
for  each  signal  is  indicated  every  30s.  Vertical  bar  represents  the  amplitude 
of  the  PAT  signal  changes  to  the  baseline. 


IV.  DISCUSSION 

Neurovascular  coupling  study  of  epilepsy  requires 
simultaneous  recording  of  neuronal  activity  and  macroscopic 
hemodynamic  signals  in  the  brain,  bringing  significant 
challenges  in  terms  of  experimental  paradigm  and  instrument 
performance.  We  successfully  integrated  the  EEG  system 
into  our  real-time  3D  PAT  system  that  is  based  on  a 
64-channel  parallel  data  acquisition  system  and  a  spherical 
transducer  array,  and  realized  the  simultaneous  recording  of 
neuron  activity  and  hemodynamic  signal  in  the  whole  rat 
brain  over  the  entire  epileptic  seizure  process. 

The  neurovascular  coupling  of  the  intraperitoneal  PTZ  rat 
model  here  is  studied  with  the  Hb02  and  HbR  signals  of  the 
SSS.  PTZ  has  been  used  experimentally  to  induce  generalized 
seizure,  where  hemodynamic  changes  can  be  found  all  over 
the  brain  [17].  As  one  of  the  main  blood  vessels  in  the  brain, 
SSS  is  highly  involved  in  the  PTZ-induced  seizures.  When 
the  seizure  begins,  we  observed  a  sudden  decrease  in  Hb02 
for  1  to  2  minutes,  followed  by  a  dramatic  increase.  The  ‘dip’ 
observed  here  is  much  similar  to  the  reported  ‘epileptic  dip’ 
with  ORIS,  but  found  on  the  cerebral  cortex.  It’s  different 
from  the  ‘initial  dip’  which  only  lasts  several  seconds.  In  the 
second  stage  of  the  seizure,  the  Hb02  signal  keeps  on 
increasing  and  finally  succeeds  the  baseline,  which  supports 
the  idea  that  the  increase  of  cerebral  blood  flow  oversupplies 
the  increased  metabolic  demands  of  the  neurons  in  epileptic 
seizures.  And  the  negative  correlation  we  observed  between 
the  Hb02  and  HbR  signals  in  SSS  here  can  be  explained  by 
the  fact  that  the  veins  don’t  usually  dilate  as  much  as  the 
arteries. 

There  are  limitations  in  our  work.  We  only  studied  the 
PAT  signals  of  SSS,  but  in-depth  information  of  the  vascular 
coupling  is  more  embedded  in  the  local  hemodynamic 
changes.  Thus  the  signals  of  different  regions  in  the  brain 


tissue  need  to  be  extracted,  based  on  which,  the  neural 
network  analysis  can  be  carried  out. 

V.  CONCLUSION 

We  for  the  first  time  present  a  combined  real-time  3D 
PAT  and  LEG  system  to  study  the  neurovascular  coupling  of 
epilepsy,  which  is  demonstrated  using  a  PTZ  induced 
generalized  seizure  model  in  rats.  PAT  data  was  collected  at 
two  different  wavelengths,  and  average  PAT  signals  of  SSS 
were  compared  with  the  corresponding  EEG  signals.  We  find 
a  long  time  ‘dip’  in  the  Hb02  signal  at  the  seizure  beginning, 
followed  by  a  profound  increase,  and  the  Hb02  and  HbR 
signals  are  negatively  correlated.  All  these  are  in  accordance 
with  other  published  results.  This  indicates  that  our  system  is 
capable  for  the  non-invasive  neurovascular  coupling  of 
epilepsy  in  small  animal  brain. 
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